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Executive  Summary 

The  research  team  made  steady  progress  in  many  areas  during  the  re¬ 
porting  period.  Review  and  research  into  the  characterization  of  the  noise 
source  continues  and  is  guiding  data  analysis  efforts.  CFD  solutions  were 
completed  and  used  to  generate  time-resolved  flow  data  for  further  analysis 
to  study  relevant  aeroacoustic  sources.  Significant  progress  was  made  on 
the  development  and  application  of  the  computational  phased  array  beam¬ 
forming  analysis.  Analysis  of  standard  stereo-PIV  data  acquired  in  October, 
2012,  was  completed  and  reduced  to  generate  a  database  for  further  study 
and  comparison  with  other  diagnostic  efforts. 
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1  Project  Objectives  and  Status 

1.1  Review  of  Program  Objectives 

The  primary  objective  of  this  research  effort  is  to  address  the  deficiency  in  the  understanding 
of  the  effect  that  near-nozzle  and  inner-nozzle  flow  conditions  have  on  jet  noise  radiation. 
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This  is  particularly  relevant  to  future  active  control  efforts  because  any  control  strategy  will 
involve  modification  of  the  boundary  conditions.  The  approach  guiding  the  current  effort  is 
three-fold: 

1.  high-fidelity  characterization  of  heated,  over-expanded  supersonic  jets, 

2.  source  identification  through  development  of  advanced  analytical  diagnostics,  and 

3.  enhanced  computational  modeling  of  hot  supersonic  jets. 

1.2  Project  Status 

An  updated  project  chart  is  shown  in  Figure  1  which  details  the  WBS  items  from  the  SOW 
and  the  completion  percentages  as  of  the  end  of  the  reporting  period.  The  research  team 
made  steady  progress  in  many  areas  during  the  reporting  period.  Review  and  research  into 
the  characterization  of  the  noise  source  continues  and  is  guiding  data  analysis  efforts.  CFD 
solutions  were  completed  and  used  to  generate  time-resolved  flow  data  for  further  analysis 
to  study  relevant  aeroacoustic  sources.  Significant  progress  was  made  on  the  development 
and  application  of  the  computational  phased  array  beamforming  analysis.  Analysis  of  stan¬ 
dard  stereo-PIV  data  acquired  in  October,  2012,  was  completed  and  reduced  to  generate 
a  database  for  further  study  and  comparison  with  other  diagnostic  efforts.  Details  of  new 
activity  and  results  are  presented  in  the  following  section. 

2  Activity  for  Current  Reporting  Period 

2.1  Significance  of  the  Jet  Conditions 

The  accurate  determination  of  the  jet  exit  conditions  is  complicated  in  general  and  more 
so  in  shock  containing  supersonic  jets.  In  the  current  study,  focused  on  conic  nozzles,  the 
non-ideal  expansion  yields  shocks  that  emanate  from  within  the  nozzle.  As  a  result,  the  exit 
plane  is  not  characterized  by  a  uniform  outflow. 

The  study  of  the  effects  of  nozzle  conditions  for  subsonic  jets  has  been  a  topic  of  discus¬ 
sion  among  a  number  of  researchers  in  the  last  decade.  In  particular,  Viswanathan[l]  did  an 
extensive  study  to  vary  Reynolds  number  and  temperature  ratio  independently  and  found  a 
transition  in  the  Reynolds  number  relevant  for  jet  noise.  It  was  concluded  that  at  the  very 
least  the  exit  conditions  are  significant  in  allowing  comparison  between  data  sets.  Likewise, 
the  more  recent  report  from  Zaman[2],  also  restricted  to  subsonic  jets,  concludes  that  vari¬ 
ations  in  the  nozzle  exit  boundary  layer  could  explain  anomalies  in  acoustic  measurements 
from  seemingly  equivalent  operating  conditions. 

For  supersonic  jets,  the  recent  report  by  Morris[3]  presented  some  comparisons  between 
model  scale  jet  noise  spectra  and  full  scale  data.  Morris  commented  that,  “The  engine  data 
does  not  roll  off  as  fast  as  the  scale  model  data  at  high  frequencies.  This  could  be  due  to  a 
number  of  reasons,  including  the  difference  in  turbulence  exiting  the  nozzle  due  to  the  engine, 
or  other  noise  sources  from  the  engine.”  This  would  suggest  that  the  correct  characterization 
of  the  nozzle  exit  conditions  is  also  of  significance  for  supersonic  jet  noise  particularly  when 
active  control  methodologies  are  being  developed  at  model  scale. 
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Figure  1.  Project  chart  showing  WBS  items  and  current  completion  status  as  of  the  end  of 
the  reporting  period. 
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Figure  2.  Location  of  the  near-field  microphones  with  comparison  to  the  locations  of  the 
Kulite  transducers  in  the  February  2012  experiment. 

A  series  of  acoustic  measurements  were  obtained  to  allow  analysis  of  some  of  these  issues. 
This  data,  acquired  in  July  2012,  includes  12  B&  K  microphones  on  a  far-field  arc  at  55  jet 
diameters  from  the  nozzle  as  well  as  10  PCB  112A22  ICP  transducers  in  the  near-field.  The 
position  of  the  PCB  microphones,  shown  in  Figure  2,  was  chosen  to  closely  follow  the  line 
of  Kulite  transducers  in  the  February  2012  experiment.  Data  was  acquired  for  two  nozzle 
configurations:  (a)  conic  nozzle  without  the  upstream  ccnterbody,  and  (b)  conic  nozzle  with 
the  upstream  centerbody.  The  stagnation  pressure  and  temperature  was  varied  to  provide  a 
range  of  NPR  from  3.9  to  6.0  and  Tjet/Too  of  2.3  to  3.5.  The  run  conditions  for  each  event 
are  provided  in  Table  1. 

The  inclusion  of  the  upstream  centerbody  in  the  nozzle  assembly  is  intended  to  replicate 
the  engine  core  in  a  military  nozzle.  The  aft  end  of  the  centerbody  is  located  approximately 
2.5  Dj  upstream  of  the  nozzle  throat.  The  wake  from  the  centerbody  yields  higher  turbulence 
levels  exiting  the  nozzle.  This  is  evident  in  the  far-field  acoustics  as  shown  in  Figure  3. 
The  increase  is  small  but  observable  in  the  peak  noise  direction  but  more  pronounced  in 
the  forward  arc.  To  further  this  analysis,  the  CFD  simulation  and  the  anticipated  PIV 
measurements  will  be  used  to  aid  in  describing  the  turbulence  levels  present  and  the  degree 
to  which  this  effects  the  noise  production  mechanisms. 

2.1.1  Near-Field,  Acoustic  Characterization 

The  near-field  acoustic  measurements  made  in  February  2012  using  the  Kulite  array  have 
been  report  previously  and  are  also  included  in  the  2012  INTERNOISE  Conference  Paper  [4], 
so  they  are  not  going  to  be  repeated  here.  The  measurements  using  the  PCB  transducers 
(Figure  2)  provide  a  second  reference  for  comparison.  Figure  4  shows  the  direct  comparison 
between  the  Kulite  and  PCB  measurements  for  the  same  jet  conditions. 

At  first  glance  the  comparison  suggests  dramatic  differences.  However,  the  measurements 
are  in  closer  agreement  when  the  influence  of  the  setup  is  considered.  The  Kulite  sensors  were 
flush  mounted  in  a  1  inch  wide  beam.  Acoustically,  the  beam  acts  as  a  finite  baffle  that  alters 
the  frequency  response  of  the  installed  sensor.  Given  the  acoustic  speed  to  be  approximately 
330  m/s,  the  width  of  the  beam  corresponds  to  a  1/2-wavelength  at  a  frequency  of  6500 
Hz.  For  a  finite  baffle,  a  pressure  doubling  is  expected  for  frequencies  higher  than  the  1/2- 
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Conic  with  Centerbody 


Date 

Patm  (PSi) 

Tj  (deg.  F) 

Po(psi) 

Mach 

NPR 

T  Ratio 

Event 

7/10/2012 

14.375562 

1398.163491 

121.042021 

1.742956 

5.214842 

3.448 

5 

7/10/2012 

14.375425 

1342.594434 

121.234251 

1.743389 

5.22159 

3.339 

6 

7/10/2012 

14.375456 

1109.583585 

121.094518 

1.741094 

5.216709 

2.901 

7 

7/10/2012 

14.374967 

1000.044245 

121.072384 

1.740227 

5.216059 

2.696 

8 

7/10/2012 

14.375403 

900.567075 

121.060892 

1.739453 

5.215546 

2.508 

9 

7/10/2012 

14.375714 

809.75283 

121.043961 

1.738722 

5.214858 

2.338 

10 

7/10/2012 

14.375367 

1204.273774 

121.012914 

1.741424 

5.213883 

3.062 

11 

7/10/2012 

14.390736 

1403.806321 

83.484092 

1.551514 

3.904342 

3.400 

12 

7/10/2012 

14.390461 

1355.352358 

83.787458 

1.552894 

3.914972 

3.310 

13 

7/10/2012 

14.389785 

1212.535094 

83.955281 

1.552588 

3.920905 

3.048 

14 

7/10/2012 

14.389661 

1098.50934 

84.197269 

1.552993 

3.929373 

2.839 

15 

7/10/2012 

14.389623 

1002.072642 

83.609451 

1.54857 

3.908896 

2.662 

16 

7/10/2012 

14.389818 

915.111415 

84.073257 

1.550565 

3.925012 

2.504 

17 

7/10/2012 

14.390414 

805.305472 

83.903022 

1.548525 

3.918973 

2.304 

18 

7/10/2012 

14.389518 

1196.742642 

85.356678 

1.560796 

3.969721 

3.016 

19 

7/10/2012 

14.383065 

1206.618868 

100.242293 

1.64306 

4.488994 

3.034 

20 

7/10/2012 

14,376516 

1194.603019 

114.391514 

1.711757 

4.983052 

3.011 

21 

7/10/2012 

14.370759 

1211.330283 

129.787898 

1.778679 

5.520823 

,  3.041 

22 

7/10/2012 

14.365041 

1211.29717 

143.121485 

1.831111 

5.987074 

3.039 

23 

Conic  without  Centerbody 


Date 

^atm  (PSi) 

T,(deg.  F) 

Po(pS') 

Mach 

NPR 

T  Ratio 

Event 

7/11/2012 

14.411241 

1406.686667 

120.721317 

1.740322 

5.193295 

3.478 

24 

7/11/2012 

14.410666 

1348.769623 

121.04801 

1.741336 

5.204776 

3.366 

25 

7/11/2012 

14.409945 

1187.632885 

121.033159 

1.740123 

5.204455 

3.061 

26 

7/11/2012 

14.410198 

1103.060283 

121 096799 

1.739783 

5.206615 

2.902 

27 

7/11/2012 

14.41005 

1013.507736 

121.332702 

1.740166 

5.214806 

2.732 

28 

7/11/2012 

14.410342 

902.982642 

121.338609 

1.739399 

5.214972 

2.525 

29 

7/11/2012 

14.409662 

802.558774 

121.694498 

1.740261 

5.227535 

2.336 

30 

7/11/2012 

14.399732 

1192.748679 

143.849986 

1.832433 

6.000418 

3.051 

31  i 

7/11/2012 

14.40591 

1203.090189 

130.160646 

1.778844 

5.52272 

3.065 

32 

7/11/2012 

14.409827 

1198.520566 

117.997351 

1.726849 

5.099088 

3.052 

33 

7/11/2012 

14.417322 

1204.078491 

100.782909 

1.644576 

4.49942 

3.057 

34 

7/11/2012 

14.423409 

1202.674135 

88.278111 

1.576668 

4.064092 

3.052 

35 

7/11/2012 

14.425184 

1396.764245 

83.84466 

1.552413 

3.909924 

3.405 

36 

7/11/2012 

14.425124 

1357.414423 

84.482442 

1.555852 

3.93206 

3.331 

37 

7/11/2012 

14.424208 

1216.224327 

84.341533 

1.553732 

3.927346 

3.072 

38 

7/11/2012 

14.424385 

1098.702075 

84.333583 

1.552593 

3.927038 

2.855 

39  i 

7/11/2012 

14.424022 

1001.747453 

84.431494 

1.5523 

3.930502 

2.678 

40 

7/11/2012 

14.423557 

798.189528 

84.067917 

1.548291 

3.918001 

2.304 

42 

7/11/2012 

14.423737 

907.192264 

84.203746 

1.550084 

3.922667 

2.504 

43 

Table  1.  Nozzle  operating  conditions  for  the  July  2012  data  set. 
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Figure  3.  Comparison  of  the  far  field  acoustics  with  and  without  the  upstream  centerbody  in 
the  nozzle  assembly. 


Figure  4.  Comparison  of  Kulite  and  PCB  measurement  of  near-field  acoustics  for  the  same 
jet  conditions 
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wavelength  equivalent  and  a  smaller  increase  for  lower  frequencies.  Looking  at  Figure  4  this 
is  exactly  the  situation.  For  frequencies  above  6500  Hz,  the  Kulite’s  report  «  6  dB  above 
the  PCB  measurement.  For  lower  frequencies,  the  Kulites  report  about  1-2  dB  above  the 
PCB.  Therefore,  the  discrepancy  is  linked  directly  to  the  frequency  response  of  the  finite 
baffle  around  the  Kulite. 

Neither  the  Kulite  nor  the  PCB  sensors  have  a  detailed  frequency  response  from  the 
manufacturer.  In  order  to  correct  the  spectral  measurements  for  the  frequency  response  of 
the  assembled  sensors,  a  separate  experiment  is  being  performed  using  a  1  inch  subsonic 
cold  jet.  The  flush  mounted  Kulite,  the  PCB  transducer,  and  a  calibration  traceable  B&K 
1/4-inch  microphone  will  be  place  as  equivalent  (x,  r )  positions  separated  by  120  degrees  in 
azimuth.  The  B&K  microphone  will  be  used  to  generate  a  frequency  response  function  for 
both  the  Kulite  and  PCB  sensors.  Finally,  the  February  2012  and  July  2012  data  will  be 
corrected  with  the  appropriate  transfer  functions  to  produce  a  final  data  set  that  will  be 
made  available  for  post  processing  and  analysis. 


2.2  Exploring  Aeroacoustic  Source  Terms 

One  of  the  inherent  features  of  an  acoustic  analogy  is  that  the  result  gives  a  representation 
of  a  non-unique  acoustic  source  that  is  consistent  with  the  measured  far-field  sound. [5] 
This  is  true  even  of  recent  computational  efforts  using  F-W&H  methods  to  propagate  near¬ 
field  pressure  fluctuations  to  a  far-field  observer  through  the  application  of  Green’s  function 
integrals.  The  drawback  is  that  much  of  the  sound  generation  is  hidden  from  exploration. 

Various  vortex  sound  methods  have  been  explored  over  the  years  as  an  alternative  or 
enhancement  to  acoustic  analogies.  The  “rational”  approach  of  Seiner  [6]  falls  into  this  cat¬ 
egory  and  suggests  a  noise  source  indicator  resulting  from  the  phase  difference  between  the 
strain  and  rotation  rates  acting  on  the  vorticity.  To  arrive  at  this  suggestion,  Seiner  applied 
a  double  divergence  directly  to  the  incompressible  Reynolds  stress  giving 


d2UiUj 

dxidxj 


-£ijk{SijU>k  +  TijUk). 


(1) 


This  neglected  variations  in  the  dynamic  pressure.  The  suggestion  Seiner  made  was  that 
any  phase  shift  between  the  strain  and  rotation  rates  in  this  sum  would  give  rise  to  sound 
production. 

To  date,  we  have  re-derived  the  expression  with  the  dynamic  pressure  included  as  pre¬ 
sented  in  the  Fievet  et  al.  abstract  for  the  2013  Aeroacoustics  conference  (included  in  the 
previous  quarterly  report): 


d2U{Uj 

dx^Xj 


=  £ijk{SijU>k  +  TijUk)  + 


(2) 


Now,  we  are  working  to  evaluate  this  expression  using  the  HRLES  CFD  data  set.  By 
developing  the  post  processing  tools  for  the  analysis,  we  are  establishing  the  building  blocks 
to  analyze  the  experimental  data  that  is  anticipated  from  the  time-resolved  PIV  experiments. 
We  are  also  establishing  the  necessary  dynamic  range  for  the  experimental  measurements  to 
allow  investigation  of  the  acoustic  source  terms. 
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Figure  5.  Illustration  of  the  internal  geometry  rendered  in  the  CFD  simulation  showing  the 
upstream  centerbody  which  replicates  the  experimental  setup. 


Conic  Md  =  1.74  Nozzle 

r„(K) 

P0  (kPa) 

Too  (K) 

Poo  (kPa) 

To/Too 

NPR 

Over-Expanded 

1005 

401.76 

300 

101.33 

3.35 

3.96 

Pressure  Matched 

1005 

531.33 

300 

101.33 

3.35 

5.24 

Ihble  2.  Stagnation  conditions  prescribed  for  the  two  CFD  simulations. 


2.2.1  Time-Resolved  CFD  Data  Set  for  Aeroacoustic  Source  Analysis 

Time  resolved  data  using  the  Hybrid  RANS/LES  CFD  simulations  have  been  compiled  for 
two  jet  conditions  consistent  with  the  experimental  measurements  ongoing  in  this  program. 
Both  simulations  use  the  conic  jet  nozzle  with  the  upstream  centerbody  as  illustrated  in 
Figure  5.  The  stagnation  conditions  prescribed  for  the  simulation  are  listed  in  Table  2. 
Additional  details  about  the  CFD  simulation  are  given  in  the  AIAA  Paper  presented  at  the 
recent  Aerospace  Sciences  Meeting.  [10] 

Having  run  the  simulation  to  convergence,  time  resolved  data  was  generated  for  post 
processing  by  saving  solution  snapshots  at  the  equivalent  of  200  kHz.  A  total  of  three 
hundred  snapshots  were  saved  for  each  jet  condition.  All  of  the  data  was  saved  in  standard 
Plot3D  format.  The  solution  files  contain  Ux,  Uy,  Uz  and  the  thermodynamic  variables  p,  P, 
T.  The  data  has  been  archived  and  made  available  to  the  other  team  members  for  further 
analysis. 

To  reduce  the  file  sizes,  the  full  CFD  solution  domain  was  not  saved  in  the  database. 
The  restricted  domain  includes  0.0  to  0.8  meters  (0  to  15.3  Dj)  in  the  axial  extent  and  out 
to  0.22  meters  (4.1  Dj)  in  the  radial  extent  as  illustrated  in  Figure  6.  This  selected  region 
represented  5  zones  from  the  CFD  solution.  To  ease  post  processing,  the  4  outer  zones  were 
combined  in  the  reduced  output  such  that  2  zones  remained:  the  inner  H-grid  in  the  jet  core 
and  an  outer  O-grid  arranged  azimuthally  as  shown  in  Figure  7 

The  acoustic  and  turbulence  terms  of  interest  include  both  mean  and  fluctuating  quan¬ 
tities.  Therefore,  a  “mean”  value  is  needed  as  the  base  flow.  The  simulation  was  not  run 
long  enough  to  generate  converged  point-statistics.  Instead  the  300  available  snapshots  were 
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Figure  6.  Illustration  of  computational  domain  saved  in  the  time-resolved  output  files  for 
further  processing. 


averaged  followed  by  azimuthal  averaging  to  generate  the  final  “mean-flow”  illustrated  in 
Figure  8. 

Nozzle  exit  velocity  profiles  were  extracted  from  the  mean  flow  calculation.  These  are  of 
interest  due  the  observation  by  many  authors  that  the  jet  noise  is  measurably  effected  by 
the  inner-nozzle  flow  conditions.  The  nozzle  exit  boundary  layer  and  turbulence  levels  are 
an  indicator  of  these  inner-nozzle  effects.  In  Figure  9  it  is  clear  that  the  presence  of  shocks 
inside  the  nozzle  disturbs  the  radial  pU  profile  in  the  over-expanded  case. 

2.2.2  Computation  of  Acoustically  Relevant  Source  Terms 

Computing  quantities  relevant  to  acoustic  generation  (i.e.  Eq.  (2))  requires  second  order 
derivatives.  A  code  for  computing  these  derivatives  was  written  based  on  the  text  of 
Anderson[ll].  This  accounts  for  the  local  divergence  of  the  grid  coordinates  in  computing 
the  derivatives  in  a  Cartesian  coordinate  system.  The  code  for  performing  the  calculations  is 
included  here  as  Appendix  A.  Subroutines  were  written  to  perform  the  necessary  calculation 
of  the  local  Jacobian  and  central  differencing.  With  these  it  is  possible  to  construct  any 
compound  derivative  of  interest.  As  an  example,  Figure  10  shows  the  result  of  computing 
£ijkSijU>k  on  one  of  the  instantaneous  snapshots.  In  the  coming  months,  this  data  will  be 
explored  further. 

2.2.3  Possible  Entropic  Effects 

Vortex  sound  based  methods  assume  incompressible  flow  which  is  not  an  applicable  assump¬ 
tion  for  the  jet  conditions  under  study  in  this  work.  Nevertheless,  it  is  still  of  value  to 
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(a)  Symmetry  Plane 


(b)  Cross-Plane  Zone  1 


(c)  Cross-Plane  Zone  2 


Figure  7.  Illustrations  of  the  computational  grid  space  retain  in  the  reduced  data  set  for 
post-processing. 


Figure  8.  Symmetry  plane  contour  of  computed  mean  axial  velocity  illustrating  the  “mean” 
flow  computed  from  the  available  300  time-resolved  snapshots. 
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Figure  9.  Nozzle  exit  profiles  of  pUx( kg/mJ  •  s)  for  the  two  NPR’s  simulated. 


Figure  10.  Snapshot  of  EijkSijUk  which  appears  to  be  largely  isolated  to  the  shear  layer. 
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evaluate  such  second-order  divergence  terms  as  they  will  still  represent  source-like  behavior. 
But  the  compressible  nature  of  the  flow  needs  to  be  better  understood.  One  possible  method 
of  addressing  this  issue  is  to  investigate  other  noise  source  effects  that  are  relevant  to  com¬ 
pressible  flows  such  as  entropy  variations.  [7]  The  entropy  contribution  has  been  demonstrated 
to  be  significant  in  some  cases.  [8,  9]  The  CFD  database  includes  the  local  thermodynamic 
variables  allowing  exploration  of  the  significance  of  these  effects  in  sound  generation  for  the 
current  jet  conditions  of  interest. 

2.3  Updates  on  Computational  Phased  Array  Work 

Significant  progress  has  been  achieved  on  the  phased  array  calculation  front  since  the  last 
quarterly  report.  The  results  were  presented  in  the  AIAA-ASM  conference  publication  and 
presentation  included  in  Appendix  C  and  D.  These  calculations  have  been  performed  on  the 
conic  nozzle  configuration  with  the  centerbody  operating  at  pressure  matched  conditions. 
One  important  result  that  emerged  from  our  phase  array  studies  was  that  it  would  be  optimal 
to  use  different  array  sizes  depending  on  the  l/3rd  octave  band  center  frequency  of  interest. 
Beamforming  calculations  at  lower  frequencies  will  make  use  of  an  array  of  larger  size  than 
calculations  at  higher  frequencies.  Since  the  CFD  calculations  have  a  high  frequency  limit 
based  on  grid  resolution  at  ~10  kHz,  we  have  selected  the  2kHz,  3.15  kHz,  5  kHz  and  8  kHz 
frequencies  as  the  focus  for  future  studies.* 

Most  beamforming  calculations  to  this  point  have  used  45  microphones  in  a  multiple-arm 
logarithmic  spiral  arrangement  (five  spiral  arms,  each  containing  nine  microphones).  How¬ 
ever,  it  was  determined  that  for  two  lower  frequencies,  (2  kHz  and  3.15  kHz),  where  a  larger 
array  size  is  warranted,  future  beamforming  calculations  will  make  use  of  99  microphones 
(nine  spiral  arms,  each  containing  11  microphones).  This  is  because  a  larger  array  requires 
more  microphones  in  order  to  avoid  aliasing  issues.  The  beamforming  calculations  at  the  5 
kHz  and  8  kHz  center  frequencies  will  use  the  45  microphone  array;  the  physical  dimension 
of  the  array  used  for  the  8  kHz  calculations  will  be  smaller  than  that  used  for  the  5  kHz 
calculation.  These  microphone  arrays  will  henceforth  be  referred  to  as  the  small,  medium 
and  large  aperture  spiral  array  (SASA,  MAS  A  and  LAS  A). 

Finally,  in  order  to  examine  the  effect  of  the  observer  viewing  angle  on  the  perceived 
location  of  the  acoustic  sources,  the  beamforming  calculations  will  be  conducted  with  phased 
arrays  rotated  to  various  azimuthal  angles.  At  each  azimuthal  location,  the  phase  center  of 
the  array  will  be  located  72  inches  from  the  aim  point;  the  aim  point  is  located  25  inches 
downstream  of  the  nozzle  exit  on  the  axis.  Figure  11  shows  the  3D  view  of  three  (out  of  9) 
of  the  azimuthally  oriented  array  layouts  (at  135°,  90°,  and  45°  to  the  upstream  jet  axis). 
The  phase  centers  for  all  9  arrays  that  will  be  considered  are  shown. 

*It  is  relevant  to  note  that  other  research  groups  performing  computational  beamforming  have  presented 
results  for  center  frequencies  as  high  as  40  kHz  without  any  verification  that  their  computational  methodology 
and  grid  are  capable  of  resolving  these  high  frequencies.  Direct  comparison  with  experimental  data  in  the 
current  work  served  as  a  validity  check  on  our  reported  results. 
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Figure  11.  Phased  arrays  rotated  at  various  angles  with  reference  to  the  jet  axis. 


2.4  Stereo  PIV  Analysis  of  Jet  Characteristics 

The  last  quarterly  report  noted  the  acquisition  of  standard  3-component  PIV  measurements 
using  the  pH  stabilized  seeding  system.  Image  data  were  acquired  using  a  pair  of  PCO.Edge 
SCMOS  cameras  with  double  exposure  capability.  The  jet  flow  was  seeded  with  0.3  micron 
alumina  particles.  Illumination  of  the  flow  was  achieved  using  a  NewWave  ND:Yag  Gemini 
PIV  laser  with  150  mJ/pulse.  Light  sheet  forming  optics  yielded  a  light  sheet  thickness  of 
approximately  1  mm.  The  ambient  air  was  not  actively  seeded;  however,  the  recirculation 
in  the  anechoic  chamber  was  sufficient  to  yield  reasonable  seed  density  in  the  entrained  air. 

The  entire  imaging  system  was  erected  on  a  precision  traverse  system  that  has  been  laser 
aligned  to  travel  parallel  to  the  geometric  axis  of  the  nozzle  system  as  shown  in  Figure  12. 
This  allowed  the  entire  measurement  entry  to  be  completed  without  any  setup  changes. 

The  experimental  measurement  involved  image  acquisition  at  four  different  axial  stations. 
Each  imaging  region  overlapped  slightly  to  allow  alignment  of  the  vector  data  during  post 
processing.  At  each  axial  station,  3  acquisition  events  were  performed  yielding  300  image 
pairs  each  to  give  a  total  of  900  image  pairs  per  station.  The  full  measurement  region  spans 
from  1.0  to  8.25  D0  in  the  axial  direction  and  out  to  0.9  Dj  in  the  radial  direction.  The 
computed  mean  velocity  for  the  entire  domain  is  shown  in  Figure  13  with  the  individual 
imaging  zones  outlined  to  show  the  overlap. 

For  this  PIV  dataset  the  conic  Mj  —  1.74  nozzle  was  used  without  the  upstream  center- 
body.  A  stagnation  temperature  of  1355  °F  and  NPR  =  3.934  were  computed  by  averaging 
the  conditions  for  all  12  image  acquisition  events.  A  standard  deviation  in  the  mean  was 
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Figure  12.  Stereo  PIV  setup  in  the  NCPA  Anechoic  Jet  Laboratory.  All  optics  are  supported 
from  the  large  breadboard  at  the  bottom  of  the  image  which  is  attached  to  a  precision  traverse 
allowing  re-positioning  with  better  than  0.01  inch  accuracy. 


Figure  13.  Mean  streamwise  velocity  obtained  from  standard  3-component  PIV  measurement 
of  the  jet  plume. 
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Figure  14.  Comparison  of  the  centerline  streamwise  velocity  between  the  HRLES  calcula¬ 
tions  and  PIV  data  for  the  conic  nozzle  without  the  centerbody  operating  at  overexpanded 
conditions. 


computed  from  the  variance  in  the  average  conditions  for  each  individual  event.  Then  an 
estimate  of  the  uncertainty  in  the  jet  conditions  was  computed  as  cr/y/N  where  N  =  12 
events.  In  this  way,  the  estimated  uncertainty  in  To  was  found  to  be  ±2  °F  and  in  NPR  was 
found  to  be  ±0.002. 

2.4.1  Preliminary  Results 

Analysis  of  the  image  data  was  completed  using  DaVis  8.  Image  pairs  were  analyzed  using 
a  multi-pass  method  with  128  x  128  pixel,  0%  overlap  square  windows  for  the  first  4  passes 
and  32  x  32  pixel,  50%  overlap  adaptive  windows  for  the  final  2  passes.  On  each  pass,  a  peak 
Q  ratio  limit  of  1.3  was  applied  and  a  two  pass  median  filter  was  used  to  remove  erroneous 
vectors.  The  final  analysis  step  applied  a  peak  Q  ratio  limit  of  1.1  and  again  applied  a  two 
pass  median  filter.  After  completing  the  image  analysis,  the  vector  data  was  post  processed 
to  center  the  vertical  extent  on  the  measured  jet  axis,  and  the  entire  data  set  was  saved  into 
Plot3D  format. 

Figure  14  shows  a  comparison  of  the  experimental  and  CFD  results  for  the  centerline 
velocity  decay.  There  is  qualitative  agreement  between  the  experiments  and  CFD,  as  far 
as  the  location  of  shock  cells  are  concerned,  for  approximately  6-8  inches  downstream  of 
the  nozzle  exit.  Farther  downstream  there  appears  to  be  something  similar  to  a  “phase- 
shifting”  between  the  simulations  and  experiments.  This  can  be  explained  by  the  progressive 
coarsening  of  the  CFD  grid  in  the  downstream  direction  and  was  expected. 

Figure  15  shows  a  comparison  of  the  mean  streamwise  velocity  profiles  at  three  axial 
stations.  The  experiments  and  simulations  show  good  agreement  with  one  another  at  the 
upstream  locations  ( x  =  3.0  inches  and  x  =  5.0  inches).  This  agreement  starts  worsening 
at  axial  stations  farther  downstream  (x  =  8.0  inches)  where  the  simulations  start  showing  a 
more  rapid  decay  of  the  streamwise  velocity  as  compared  to  the  experiments  (again  due  to 
the  grid  coarsening).  The  velocity  defect  at  the  centerline,  that  is  seen  in  both  experiments 
as  well  as  the  simulations,  is  a  result  of  the  slip  plane. 
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Figure  15.  Comparison  of  the  mean  streamwise  velocity  profile  at  various  axial  locations 
between  the  HRLES  calculations  and  PIV  data  for  the  conic  nozzle  without  the  centerbody 
operating  at  overexpanded  conditions. 

2.5  Progress  on  HDR-PIV  Diagnostic  Efforts 

2.5.1  DEVOLS  Processing  for  HDR-PIV 

A  method  for  obtaining  enhanced  signal-to-noise  ratios  with  the  high-dynamic-range  particle 
image  velocimetry  (HDR-PIV)  has  been  established.  This  dynamic  evaluation  via  ordinary 
least  squares  (DEVOLS)  method  has  demonstrated  measured  improvement  over  standard 
PIV  for  a  simulated  data  set.  Figure  16  shows  a  measurable  improvement  over  standard  PIV 
analysis  for  a  simulated  Hamel-Oseen  vortex.  Complete  details  of  the  DEVOLS  method  and 
results  were  presented  at  the  recent  2013  AIAA  Aerospace  Sciences  Meeting  (Paper  2013- 
0774  attached  in  Appendix  E). 

2.5.2  Current  Status  of  Experimental  Efforts 

Efforts  in  Year  1  to  acquire  time-resolve  PIV  data  were  unsuccessful  due  to  unanticipated 
operational  issues  with  the  Cordin  camera  system.  The  camera  was  sent  to  Cordin  for  service 
and  was  returned  to  Auburn  U.  in  December,  2012.  The  pulse-burst  laser  system  and  camera 
were  brought  to  NCPA  by  Auburn  graduate  student  Bryan  Brock  to  initiate  Year  2  testing 
during  the  Spring  months.  Setup  of  the  experiment  is  currently  underway. 

The  time-resolve  CFD  data  was  can  serve  as  a  means  to  develop  analysis  methods  that 
will  be  used  on  the  high-frame-rate  PIV  data  anticipated  to  be  generated  in  the  coming 
months.  As  such,  the  CFD  solution  file  was  used  to  select  a  region  of  interest  near  the  nozzle 
exit  and  generate  simulated  PIV  data. 

The  measurement  domain  shown  in  Figure  17  was  chosen  as  a  trade-off  between  imaging 
capability  and  measurement  goals.  As  noted  in  Section  1.1,  the  near-nozzle  flow  conditions 
are  of  primary  interest  in  the  current  effort.  To  maximize  the  optical  setup,  a  1.5  inch  field  of 
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(a)  Standard  PIV  Processing  (b)  DEVOLS  PIV  Processing 

Figure  16.  Illustration  of  the  reduced  absolute  error  afforded  by  application  of  DEVOLS 
processing  to  a  simulated  time-resolved  PIV  data  set. 


view  was  chosen  giving  a  suitable  spatial  resolution  to  resolve  the  high-shear  layer  gradients 
in  this  region  of  the  jet.  Because  of  the  high-shear  rates  in  the  region  of  interest,  the  high- 
dynamic-range  analysis  methods  developed  in  Year  1  will  be  well  suited  for  analyzing  the 
image  data,  and  the  analysis  of  any  phase  misalignment  between  the  shear  and  rotation  rate 
components  of  the  noise  source  should  be  observable. 

3  Technical/Cost  Status  &  Problem  Areas 

3.1  Status  Summary 

Overall,  the  status  of  the  program  is  on  track  and  progress  is  being  made  in  all  areas  of  the 
research  effort.  As  we  generate  experimental  data  and  develop  the  diagnostic  efforts,  oppor¬ 
tunities  for  collaboration  with  other  jet  noise  reduction  research  teams  are  being  identified. 
We  have  had  discussions  with  Dr.  Kent  Gee  concerning  the  sharing  of  our  acoustic  data  for 
comparison  to  his  full  scale  acoustic  data.  We  have  also  had  discussions  with  Dr.  Todd  Lowe 
concerning  the  possibility  of  utilizing  his  laser  diagnostic  system  to  measure  the  jet  exhaust 
characteristics  in  the  NCPA  Anechoic  Jet  Lab.  These  collaborative  opportunities  will  be 
developed  further  during  the  next  reporting  period. 

3.2  High- Frame-Rate  Camera  for  HDR-PIV 

Following  service  by  Cordin,  the  camera  was  brought  back  to  NCPA  in  January  2013  for  a 
second  attempt  at  data  collection.  Initial  results  with  the  camera  are  positive.  The  software 
issues  that  previous  caused  the  system  to  lock-up  have  been  eliminated,  and  the  relative  focus 
of  the  8  individual  imaging  units  seems  to  be  in  closer  agreement.  The  aperture  limitations 
of  the  optical  system  are  still  an  inherent  drawback  of  the  system.  A  large  aperature  135mm 
lens  has  been  purchased  and  is  currently  being  used  as  the  primary  lens  allowing  greater 
light  collection.  Currently  the  system  is  being  readied  for  measurement  in  the  next  week. 
By  performing  measurements  early  in  the  semester  we  hope  to  identify  any  remaining  issues 
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(a)  Region  of  Interest 


(b)  Simulated  PIV  Data 


Figure  17.  Illustration  of  the  simulated  PIV  data  extracted  from  the  CFD  solution. 


and  be  able  to  address  them  in  order  to  guarantee  a  successful  synchronized  PIV / Acoustics 
measurement  later  in  the  Spring. 

Contact  has  been  made  with  the  manufacturer  of  a  new  high  speed  camera,  the  Kirana. 
Its  capabilities  exceed  that  of  the  Cordin.  The  manufacturer  has  responded  with  great 
interest  in  supporting  a  demonstration  of  the  Kirana  system  in  April  which  will  provide 
us  with  an  opportunity  to  acquire  time-resolved  PIV  measurements  with  a  second  imaging 
system. 

3.3  Other  Potential  Imaging  Diagnostic  Options 

Plenoptic  imaging  for  volumetric  PIV  is  a  technology  that  is  in  development  by  our  Auburn 
partners.  Measurements  were  obtained  with  Auburn’s  Plenoptic  camera  last  Fall.  Initial 
results  of  the  analysis  are  very  promising.  As  a  result,  we  axe  planning  to  perform  a  second 
experiment  with  the  Plenoptic  PIV  system  during  the  next  reporting  period  with  more 
thorough  setup  and  calibration  allowing  a  detailed  analysis  of  the  flow. 

4  Publications,  Meetings,  and/or  Travel 

4.1  Running  List  of  Publications  Produced 

•  Baars,  W.,  Tinney,  C.,  and  Wochner,  M.  “Nonlinear  Noise  Propagation  from  a  Fully 
Expanded  Mach  3  Jet.”  50th  AIAA  Aerospace  Sciences  Meeting ,  Paper  2012-1177. 

•  Murray,  N.,  Lyons,  G.,  Tinney,  C.,  Donald,  B.,  Baars,  W.,  Thurow,  B.,  Haynes,  H.,  and 
Panickar,  P.  “A  Laboratory  Framework  for  Synchronous  Near/Far-Field  Acoustics  and 
MHz  PIV  in  High-Temperature,  Shock- Containing  Jets.”  Proceedings  of  the  Intemoise 
2012/ASME  NCAD  Meeting ,  Paper  ASME/NCAD-1270. 
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•  Baars,  W.  and  Tinney,  C.  “Scaling  Model  for  Nonlinear  Supersonic  Jet  Noise.”  Bulletin 
of  the  American  Physical  Society ,  57:17,  Abstract  D24:8,  San  Diego,  CA.  http:// 
meetings .  aps .  org/link/BAPS .  2012 .  DFD . D24 . 8.  [included  in  Appendix  B] 

•  Panickar,  P.,  Erwin,  J.,  Sinha,  N.,  Murray,  N,  and  Lyons,  G.  “Localization  of  Acoustic 
Sources  in  Shock-Containing  Jet  Flows  Using  Phased  Array  Measurements.”  51st 
AIAA  Aerospace  Sciences  Meeting,  Paper  2013-0613.  [included  in  Appendix  C  and  D] 

•  Haynes,  R.,  Brock,  B.,  and  Thurow,  B.  “Application  of  MHz  Frame  Rate,  High 
Dynamic  Range  PIV  to  a  High-Temperature,  Shock-Containing  Jet.”  51st  AIAA 
Aerospace  Sciences  Meeting ,  Paper  2013-0774.  [included  in  appendix  E] 

•  Fievet,  R.,  Tinney,  C.,  Murray,  N.,  Lyons,  G.,  and  Panickar,  P.  “Acoustic  Source 
Indicators  using  LES  in  a  Fully  Expanded  and  Heated  Supersonic  Jet.”  accepted  for 
presentation  at  the  2013  AIAA  Aeroacoustics  Conference. 

•  Baars,  W.,  and  Tinney,  C.  “Quantifying  Shock-Type  Acoustic  Waveform  Structures 
Emitted  by  a  Fully-Expanded  Mach  3  Jet.”  accepted  for  presentation  at  the  2013 
AIAA  Aeroacoustics  Conference. 

•  Baars,  W.  and  Tinney,  C.  “Quantifying  Supersonic  Jet  Noise  Crackle.”  accepted  for 
presentation  at  the  2nd  Symposium  on  Fluid- Structure- Sound  Interactions  and  Control, 
2013. 

•  Baars,  W.,  Tinney,  C.,  and  Wochner,  M.  “Nonlinear  Distortion  of  Acoustic  Waveforms 
from  High-Speed  Jets.”  J.  Fluid  Mech.  [in  review] 

5  Planned  Activities  for  Next  Reporting  Period 

•  Continue  analysis  of  noise  source  indicators  using  the  computational  simulation  data. 

•  Bring  the  uncertainty  analysis  efforts  up  to  date. 

•  Conduct  year  2  experimental  efforts  to  measure  acoustic  source  quantities. 

•  Reach  out  to  other  JNR  program  participants  to  foster  potential  collaboration. 
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Appendix  A:  FORTRAN  Code  for  Computing  Derivatives  from  CFD  Data 
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computingDerivatives_l_001 . f  90 
Version  1.001 

I  have  attempted  to  generalize  the  process  of  computing  derivatives  on  a 
computational  grid  space.  This  requires  the  calculation  of  the  Jacobian 
of  the  transformation  between  grid  units  and  physical  coordinates.  By 
placing  all  the  specifics  in  subroutines  and  functions,  the  effort  of 
calculating  a  derivate  is  reduced  to  specifying  the  variable  of  interest 
and  the  coordinate  direction  in  which  the  partial  derivative  is  desired. 
Now,  the  function  calls  can  be  used  to  construct  calculations  involving 
multiple  derivatives. 

Thanks  to  James  Erwin  (CRAFT  Tech)  for  giving  me  the  building  blocks 
to  build  this  code. 

Details  of  the  computation  of  the  grid  metrics  using  the  Jacobian  can 
be  found  in  the  CFD  text  by  Anderson. 

The  remaining  computation  of  interest  would  be  second  partials.  This 

code  can  easily  be  used  to  construct  mixed  second-partials 

(i.e.  dA2/dxdy) ,  but  a  savings  in  computational  noise  may  be  gained 

if  second-partials  such  as  dA2/dxA2  were  computed  in  a  single 

step.  However,  this  would  require  the  addition  of  a  new  function  that 

would  (a)  compute  derivatives  of  the  local  grid  metrics  AND  (b)  perform  the 

1st  order  2nd-central-dif ference  .  .  . 

i.e.  (  f (i+1 )  +  f(i-l)  -  2f(i)  J  /  [  dx~2  ] 

This  step  has  not  been  completed  in  this  version.  For  now,  second-partials 
can  be  obtained  by  twice  applying  the  partialOne(]  function. 

(28  Jan.  2013  ::  Nathan  Murray) 


PROGRAM  computingDerivatives 
IMPLICIT  NONE 
TYPE  func 

REAL, DIMENSION (:,:,:,:),  POINTER  ::  f  lean  be  either  grid  or  function  data 
INTEGER: :ni,nj,  nk,nf 
END  TYPE 

TYPE  neighborlDs 

INTEGER  : :  zID  !  Zone  number  of  neighbor 
INTEGER  : :  ii  !  i-index  of  neighbor 

INTEGER  ::  jj  !  j-index  of  neighbor 

INTEGER  ; :  kk  !  k-index  of  neighbor 

LOGICAL  isBoundary  !  TRUE  if  neighbor  is  the  boundary 
!  FALSE  for  internal  nodes 

END  TYPE 

INTEGER  ::  i, j , k, n, ng, ni, n j , nk, nvar, f 
INTEGER  ::  ii,jj,kk 
INTEGER  ::  ll,mm,nn 

INTEGER, DIMENSION ( 3 ) , PARAMETER  ::  UVWlD«  (/2, 3,4/) 

TYPE (func) , DIMENSION ( : )  ,  POINTER  ::  grO  !  the  grid 

TYPE (func) , DIMENSION (:),  POINTER  ::  gfun  !  the  flow  variables 

TYPE (func) , DIMENSION (:), POINTER  ::  dgfun  !  the  computed  derivative 

TYPE (neighborlDs) , DIMENSION (3,2)  : :  nID 
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64 

65 

66 

67 

68 

LOGICAL  : :  onBoundary 

CHARACTER (1024)  ::  fn 

WRITE (* , 1 (A) ' )  "Reading  Data" 

69 

70 

WRITE (fn, ' (A) ' )  "VolumeData  ConicWithCB  OverExp/SubsetXYZFile .bin" 

71 

OPEN (UNIT-101, FILE-TRIM (fn) , STATUS-"old" , FORM-" unformatted" , ACTION-" read" ) 

72 

READ (101)  ng 

73 

WRITE ( *, ' (A, 112) 1 )  "Number  of  input  zones  =",  ng 

74 

ALLOCATE (grO (ng) ) 

75 

READ (10 1 )  (grO (n) %ni, grO (n) %n j , grO (n) %nk, n-1, ng) 

76 

DO  n*l,ng 

77 

ALLOCATE (grO (n) %f (grO (n)  %ni , grO (n) %nj , grO  (n) %nk, 3) ) 

78 

READ (101)  grO  (n) %f 

79 

ENDDO 

80 

CLOSE ( 101 ) 

81 

WRITE (*,  ' (A)  * )  "Grid  File  Loaded" 

82 

83 

WRITE (f n, * (A) 1 )  "VolumeData  ConicWithCB  OverExp/Snapshot  001. bin" 

84 

OPEN (unit=101 , f ile=TRIM (fn) , FORM= ' unformatted 1 , STATUS*' old  * , ACTION=  *  read' ) 

85 

READ (101 )  ng 

86 

IF  (  ng  .NE.  SIZE(grO))  STOP  "grid  and  function  file  not  same  number  of  blocks" 

87 

ALLOCATE (gfun (ng)  ) 

88 

READ (101)  (gfun(n) %ni,gfun(n) %nj,gfun(n) %nk,gfun(n) %nf,n-l,ng) 

89 

DO  n=l,ng 

90 

ALLOCATE (gfun (n) %f (gfun (n) %ni, gfun (n) %nj , gfun (n) %nk,gfun (n) %nf ) ) 

91 

READ ( 101 )  gfun (n) %f 

92 

END  DO 

93 

|  CLOSE (101) 

94 

WRITE ( *,  ' (A) ' )  "Solution  File  Loaded" 

95 

96 

97 

98 

99 
100 

J  Convert  the  (x,y,z)  grid  data  to  (x,r,phi)  ...  READ  discussion  of  issues 

J  in  the  subroutine  below  before  using  this  conversion  . . . 
i !  CALL  convertXYZtoXRPhi (grO) 

!  Allocate  array  for  output  of  calculation: 

101 

1  !  Variable  1  ::  du/dz 

102 

?  Variable  2  ::  epsilon_[i jk] (s_[i j ] *omega_ [k] ) 

103 

ALLOCATE (dgfun (ng) ) 

104 

|  DO  n*l,ng 

105 

dgfun (n) %ni=gfun (n) %ni 

106 

dgfun (n) %nj=gfun (n) %nj 

107 

dgfun  (n) %nk=gfun (n) %nk 

108 

dgfun (n) %nf-2 

109 

ALLOCATE (dgfun (n) %f (dgfun (n) %ni, dgfun (n) %nj , dgfun (n) %nk,  dgfun  (n) %nf ) ) 

110 

ENDDO 

111 

112 

DO  n-l,ng 

113 

WRITE ( * , ’ (A,  15)  '  )  "Computing  Zone  ",n 

114 

DO  i=l,gr0 (n) %ni 

115 

DO  j=l , grO (n) %nj 

116 

DO  k=l, grO (n) %nk 

117 

J  Find  indices  for  my  neighbors  ->  nID 

118 

CALL  setNeighbors (n, i, j , k,gr0  (n) %ni,gr0 (n) %nj,gr0 (n)  %nk,nID) 

119 

120 

!  Ask  for  derivative  of  variable  in  p  coordinate  direction  . . . 

121 

!  In  the  current  dataset  the  streamwise  velocity  is  variable  2  ... 

122 

J  partialOne (dataSet, varlD, gridData, i, j , k, neighborlDs , derivDirection) 

123 

dgfun (n) %f (i, j , k, 1)  =  partialOne (gfun, 2, gr0,i, j,  k,nID, 3) 

124 

J  first  partial  of  variable  2  in  gfun  w.r.t.  direction  3  in  grO  ... 

125 

!  this  is  du  /  dz 

126 
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127 

!  Construct  the  omega  dot  strain-rate  tensor  product  similar  to 

128 

!  Seiner's  rational  approach  ...  however,  this  is  NOT  the  Seiner 

129 

!  term  of  interest  ...  here  I'm  computing  the  TOTAL  not  the 

130 

!  FLUCTUATING  quantity  that  is  really  of  interest. 

131 

dgfun (n) %f (i , j , k, 2)  -  0.0 

132 

!  Loop  over  (ll,mm,nn)  *  1  ->  3  to  perform  tensor  product 

133 

DO  11-1,3;  DO  mm-1,3;  DO  nn«l,3 

134 

IF  (  pt (11, mm, nn)  .EQ.  0.0  )  THEN 

135 

!  pt  —  0.0  is  the  permutation  tensor  =  0  so  we  can  cycle 

136 

!  to  the  next  interation  of  the  DO  loop  ... 

137 

CYCLE 

138 

ELSE 

139 

dgfun (n) %f (i, j , k, 2)  -  dgfun (n) %f (i, j  ,  k,  2)  +  & 

140 

pt(ll,mm,nn)  *  0.5  *  omega (gfun, uvwID, grO,  i,  j  ,  k,  nID, nn)  *  & 

141 

(  partialOne (gfun, uvwID (11) , grO, i, j , k, nID, mm)  +  & 

142 

partialOne (gfun, uvwID  (mm) , grO, i, j , k, nID,  11)  ) 

143 

!  The  above  sums  up  all  the  terms  of 

144 

!  epsilon_lmn  *  s_lm  *  omega_n 

145 

)  where  s  lm  is  computed  as 

146 

!  0.5  *  (  du  1/dx  m  +  du  m/dx  1  ) 

147 

ENDIF 

148 

ENDDO;  ENDDO;  ENDDO 

149 

150 

ENDDO 

151 

ENDDO 

152 

ENDDO 

153 

ENDDO 

154 

155 

WRITE (*, ' (A) ' )  "Writing  Output" 

156 

157 

WRITE (fn, ' (A) ')  "DerivativesOutput.bin" 

158 

OPEN (UNIT-101,  FILE-TRIM (fn) , FORM- 1  unformatted' , STATUS- ’ replace ' , ACTION- ' write ' ) 

159 

WRITE (101)  ng 

160 

WRITE ( 101 )  (dgfun (n) %ni, dgfun (n) %n j , dgfun  (n) %nk, dgfun (n) %nf ,  n«l , ng) 

161 

DO  n«l,ng 

162 

WRITE (101)  dgfun (n) %f 

163 

ENDDO 

164 

CLOSE ( 101 ) 

165 

166 

WRITE (*, ' (A) 1 )  "Done" 

167 

168 

j  CONTAINS 

169 

170 

SUBROUTINE  convertXYZtoXRPhi (grid) 

171 

!  convertXYZtoXRPhi  assumes  'grid'  contains  a  cartesian  grid  with  variables 

172 

!  ,x,y,z'  and  converts  to  a  cylindrical  grid  with  'x^/phi' 

173 

! -  Discussion  - 

174 

!  This  conversion  assumes  a  right-handed  coordinate  system  where  'x'  is 

175 

I  the  axis  of  the  cylindrical  coordinate  system.  The  calculation  of  1  r* 

176 

!  is  then  straightforward  as  (yA2  +  zA2)A0.5.  However,  the  determination 

177 

!  of  the  azimuthal  angle  phi  requires  care.  The  FORTRAN  intrinsic  ATAN2  is 

178 

!  aware  of  this  issue  and  takes  into  account  what  quatrant  (y,z)  is  in 

179 

!  when  returning  the  angle.  However,  the  result  is  restricted  to  the 

180 

!  range  -pi/2  ->  +pi/2  where  the  angle  is  measured  as  the  counter-clockwise 

181 

i  !  rotation  from  the  positive  1 y '  axis.  NOTICE  —  If  this  conversion  is 

182 

!  performed  prior  to  computing  derivatives  there  WILL  BE  ambiguity  at  the 

183 

J  negative  'y'  axis  where  a  "jump"  will  appear  from  +pi/2  to  -pi/2  between 

184 

!  neighboring  grid  points. 

185 

j -  Concerning  the  Present  Situation  - 

186 

)  I  have  created  this  problem  in  some  degree  by  generating  'neighbor'  rules 

187 

!  that  allow  the  azimuthal  grid  index  to  wrap  . . .  one  way  of  avoiding  the 

188 

L  !  above  problem  is  to  compute  only  forward  or  backward  derivatives  at  the 

189 

!  negative  * y '  axis  ...  in  essence  setting  the  negative  'y'  axis  as  an 
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190 

!  actual  boundary.  Another  possible  solution  would  be  the  assignment  of 

191 

!  a  separate  type  of  boundary  ...  a  periodic  boundary  at  the  negative  *  y' 

192 

!  axis  . . .  but  this  would  require  a  separate  set  of  conditions  to  handle 

193 

!  central-differences  across  the  periodic  boundary  ...  for  now,  I'm  not  going 

194 

!  to  implement  either  of  these  solutions  but  only  mention  them  here  for 

195 

!  reference. 

196 

TYPE ( func ) , DIMENSION ( : ) , INTENT ( INOUT)  ::  grid 

197 

INTEGER  ::  i,  j  ,  k 

198 

REAL  : :  r,phi, y, z 

199 

200 

DO  n*l , SIZE (grid) 

201 

DO  i*l,grid (n) %ni 

202 

DO  j=l,grid (n) %nj 

203 

DO  k-1 , grid (n) %nk 

204 

y  *  grid (n) %f (i, j ,k, 2) ; z  «  grid (n) %f (i, j , k, 3) 

205 

grid (n) %f (i, j , k, 2 )  *  SQRT(y**2.0  +  z**2.0) 

206 

grid (n) %f (i, j  ,  k,  3 )  «  ATAN2 ( z, y) 

207 

ENDDO 

208 

ENDDO 

209 

ENDDO 

210 

ENDDO 

211 

212 

END  SUBROUTINE  convertXYZtoXRPhi 

213 

214 

SUBROUTINE  setNeighbors (currentZone, ii, j  j , kk, imax, jmax, kmax, neighborlD) 

215 

!  setNeighbors  builds  the  small  array  'neighborlD'  which  identifies  the  grid 

216 

!  location  of  all  6  neighbors  to  the  current  cell.  These  are  subsequently  used 

217 

•  !  to  compute  derivatives  as  desired.  There  are  some  calculations  here  that  are 

218 

!  strictly  dependent  on  the  current  data  set.  The  idea  was  to  use  this 

219 

!  subroutine  to  identify  the  neighbors  so  that  the  computation  of  derivatives 

220 

!  could  be  completely  general . 

221 

INTEGER, INTENT (IN)  ::  currentZone, ii, j j , kk, imax, jmax, kmax 

222 

TYPE (neighbor IDs), DIMENSION INTENT  (OUT)  ::  neighborlD 

223 

!  neighborlD  [  [i,j,k  Direction],  (behind, ahead]  ] 

224 

225 

!  Initialize  all  'isBoundary'  logicals  to  FALSE  to  ensure  they  are  only 

226 

I  !  set  to  true  if  necessary. 

227 

neighborlD ( : , : ) %isBoundary  =  .FALSE. 

228 

229 

230 

!  Set  neighbors  in  the  i-Direction 

231 

232 

!  The  ’i’  direction  is  the  axial  direction  and  is  continuous 

233 

!  without  any  zonal  boundaries  in  the  current  data  set.  All 

234 

!  I  have  to  check  for  is  the  grid  space  boundaries. 

235 

!  I  can  do  this  using  the  MIN  and  MAX  commands  without 

236 

! requiring  IF  THEN  blocks. 

237 

neighborlD (1, 1 : 2) %zID  =  currentZone 

238 

neighborlD (1, 1  ) %ii  =  MAX(ii-l,l) 

239 

neighborID(l,  2  )  %ii  *=  MIN  (ii+1 ,  imax) 

240 

neighborlD ( 1, 1 : 2 ) %j j  «  jj 

241 

neighborlD (1, 1 : 2) %kk  «  kk 

242 

IF  (  ii  .EQ.  1  )  THEN 

243 

neighborlD (1, 1 ) %isBoundary  *  .TRUE. 

244 

ELSE  IF  (  ii  .EQ.  imax  )  THEN 

245 

neighborlD (1, 2 ) %isBoundary  =  .TRUE. 

246 

ELSE 

247 

neighborlD (1 , 2 ) %isBoundary  =  .FALSE. 

248 

ENDIF 

249 

250 

!  In  the  ’j’  and  'k1  directions,  all  my  neighbors  have  the 

251 

!  same  ii  index,  so  I  can  set  that  here  . . . 

252 

neighborID(2 : 3, 1:2) %ii  =  ii 
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253 

254 

IF  (  currentZone  .EQ.  1  )  THEN 

255 

256 

!  Set  Neighbors  in  Zone  1  ...  the  outer  zone  . . . 

257 

!  The  ' j '  direction  is  the  radial  direction  ... 

258 

!  The  ' k  *  direction  is  the  azimuthal  direction  ... 

259 

i _ 

2  60 

!  The  'ahead*  j-neighbor  is  always  in  the  current  zone 

261 

neighbor ID (2, 2) %zID  -  currentZone 

262 

neighbor ID (2, 2)  %kk  -  kk 

263 

!  ...  but  may  be  on  the  outer  boundary 

264 

neighbor ID (2, 2) % j  j  -  MIN ( j j  +  1 , jmax) 

265 

IF  (  jj  .EQ.  jmax  )  neighborlD (2, 2 ) %isBoundary  -  .TRUE. 

266 

267 

!  If  I'm  not  on  the  inner  boundary  ... 

268 

IF  (  jj  .GT.  1  )  THEN 

269 

!  the  'behind'  neighbor  is  straightforward 

270 

neighborlD (2, 1 ) %zID  *  currentZone 

271 

neighbor ID (2, 1 ) % j  j  -  jj-1 

272 

neighborlD (2, 1) %kk  -  kk 

273 

ELSE 

274 

!  The  'behind'  neighbor  is  in  zone  2  ...  the  inner  zone 

275 

neighborlD (2, 1 ) %zID  *  2 

276 

!  The  specific  neighbor  now  depends  on  which  sector  I 

277 

!  am  in,  A  simple  INTO  function  can  indicate  the  sector 

278 

!  given  that  there  are  30  elements  in  each  of  the  4  sectors. 

279 

IF  (  INT  (  ( kk- 1 )  /  30 )  +1  .EQ.  1  )  THEN  .'  sector  1 

280 

neighborlD (2, 1 ) % j  j  -  1 

281 

neighborlD (2, 1 ) %kk  -  30  -  kk  +  1 

282 

ELSE  IF  (  INT  (  (kk-1) /30) +1  .EQ.  2  )  THEN  !  sector  2 

283 

neighborlD (2, 1 ) % j j  -  kk  -  30 

284 

neighborlD (2, 1 ) %kk  -  1 

285 

ELSE  IF  (  INT ( (kk-l)/30)+l  . EQ .  3  )  THEN  !  sector  3 

286 

neighborlD (2, 1 ) % j  j  -  30 

287 

neighborlD (2, 1) %kk  -  kk  -  60 

288 

ELSE  IF  (  INT  ( (kk-l ) /30 ) +1  .EQ.  4  )  THEN  !  sector  4 

289 

neighborlD (2 , 1) % j j  -  30  -  (kk  -  90)  +  1 

290 

neighborlD (2 , 1 )  %kk  -  30 

291 

ENDIF 

292 

ENDIF 

293 

294 

!  Set  the  neighbors  in  the  * k*  direction  ... 

295 

!  Both  'ahead*  and  'behind'  neighbors  are  in  the  same  zone  ... 

296 

neighborlD  (3, 1 : 2) %zID  -  currentZone 

297 

neighborlD (3, 1 : 2) %j j  *  jj 

298 

J  ...  but  there  is  a  wrap-around  in  the  azimuthal  direction 

299 

neighborlD  (3, 1) %kk  -  kk-1 

300 

neighborlD (3, 2) %kk  -  kk+1 

301 

IF  (  kk  .EQ.  1  )  THEN 

302 

neighborlD (3, 1 ) %kk  «  kmax  !  'behing'  neighbor  has  to  wrap 

303 

ELSE  IF  (  kk  .EQ.  kmax  )  THEN 

304 

neighborlD (3 , 2) %kk  -  1  !  'ahead'  neighbor  has  to  wrap 

305 

ENDIF 

306 

ELSE  IF  (  currentZone  .EQ.  2  )  THEN 

307 

308 

!  Set  Neighbors  in  Zone  2  ...  the  inner  zone  . . . 

309 

!  'j'  and  'k'  represent  a  curvilinear  coordinate  system 

310 

311 

!  In  zone  2  I  have  to  be  careful  with  the  boundaries 

312 

!  Set  the  j-Direction  neighbors  first 

313 

IF  (  jj  .EQ.  1  )  THEN 

314 

neighborlD (2, 1 ) %zID  -  1  !  behind  neighbor  in  Zone  1 

315 

neighborlD (2, 1 ) %j j  -  1  )  All  Zone  1  neighbors  to  Zone  2  have  j«l 
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316 

neighborlD (2, 1) %kk  »  30  -  kk  +  1 

317 

neighborlD (2, 2 ) %zID  »  currentZone 

318 

neighbor ID (2, 2 ) % j  j  =  jj  +  1 

319 

neighbor ID (2/  2)  %kk  -  kk 

320 

ELSE  IF  (  jj  .EQ.  jmax  )  THEN 

321 

neighborlD (2, 1) %zlD  =  currentZone 

322 

neighbor ID (2, 1 ) % j j  =  jj-1 

323 

neighborlD  (2, 1)  %kk  *=  kk 

324 

neighborlD (2, 2) %zID  »  1  !  ahead  neighbor  in  Zone  1 

325 

neighborlD (2, 2) %jj  -  1  !  All  Zone  1  neighbors  to  Zone  2  have  j=l 

326 

neighborlD (2, 2) %kk  »  kk  +  60 

327 

ELSE  !  I'm  not  on  a  boundary  and  the  neighbors  are  straightforward 

328 

neighborlD (2, 1 : 2 ) %zID  *  currentZone 

329 

neighborlD (2 , 1  )%jj  =  jj-1 

330 

neighborlD (2, 2  ) % j  j  =  jj  +  1 

331 

neighborlD (2,  1 :2) %kk  »  kk 

332 

ENDIF 

333 

334 

!  Set  the  k-Direction  neighbors 

335 

IF  (  kk  .EQ.  1  )  THEN 

336 

neighborlD (3/ 1)  %zID  *=  1  !  behind  neighbor  in  Zone  1 

337 

neighborlD (3, 1) %jj  =  1 

338 

neighborlD (3, 1 ) %kk  =  30  +  jj 

339 

neighborlD (3/ 2) %zID  =  currentZone 

340 

neighborlD  (3,  2)  %jj  *=  jj 

341 

neighborlD  (3,  2)  %kk  *=  kk+1 

342 

ELSE  IF  (  kk  .EQ.  kmax)  THEN 

343 

neighborlD (3, 1) %zID  =  currentZone 

344 

neighborlD (3, 1 ) % j  j  -  jj 

345 

neighborlD (3, 1 ) %kk  =  kk-1 

346 

neighborlD (3, 2) %zID  «  1  !  ahead  neighbor  in  Zone  1 

347 

neighborlD  (3, 2 )  %  j  j  *=  1 

348 

neighborlD (3, 2 ) %kk  *=  90  +  30  -  jj  +  1 

349 

ELSE  !  I'm  not  on  a  boundary  and  the  neighbors  are  straightforward 

350 

neighborID(3, 1 :2) %zID  *=  currentZone 

351 

neighborlD (3, 1 : 2 ) %j j  =  jj 

352 

neighborlD (3, 1  )%kk  -  kk-1 

353 

neighborlD (3, 2  ) %kk  -  kk+1 

354 

ENDIF 

355 

ELSE 

356 

ENDIF 

357 

358 

END  SUBROUTINE  setNeighbors 

359 

360 

FUNCTION  partialOne (vFunc, nVar, coord, i, j , k, nID, dir)  RESULT (varDer) 

361 

!  partialOne  computes  the  first  derivative  of  variable  'nVar*  in  'vFunc' 

362 

!  at  grid  location  (i,j,k)  in  the  'dir'  coordinate  direction 

363 

!  using  the  neighbors  defined  in  nID.  The  coordinate  direction 

364 

!  corresponds  to  one  of  the  grid  axes  ...  i.e.  (x,y,z)  or  (x,r,phi) 

365 

IMPLICIT  NONE 

366 

TYPE  (func) , DIMENSION ( : ) , INTENT ( IN)  ::  vFunc 

367 

TYPE (func), DIMENSION (:), INTENT (IN)  ::  coord 

368 

INTEGER, INTENT (IN)  ::  i , j , k, dir , nVar 

369 

TYPE (neighborlDs) , DIMENSION  (:,:), INTENT (IN)  ::  nID 

370 

REAL  ::  varDer, dl , d2 , d3 

371 

372 

REAL, DIMENSION (3,  3)  ::  gMets 

373 

374 

gMets  =  setMetrics (coord, i, j , k, nID) 

375 

376 

!  Differencing  of  the  function  variable  on  the  computational  space 

377 

dl  -  dol (vFunc, nVar, nID, 1) 

378 

d2  =  dol (vFunc, nVar, nID, 2) 
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379 

d3  -  dol (vFunc, nVar, nID, 3) 

380 

381 

varDer  -  gMets (1, dir ) *dl  +  gMets (2, dir ) *d2  +  gMets (3, dir ) *d3 

382 

383 

END  FUNCTION  partialOne 

384 

385 

FUNCTION  setMetrics (coord, ii, j j , kk, nID)  RESULT (gridMet) 

386 

!  setMetrics  returns  the  3x3  matrix  of  grid  metrics  necessary  for 

387 

!  computing  derivatives. 

388 

IMPLICIT  NONE 

389 

TYPE (func) , DIMENSION ( : )  , INTENT (IN)  ::  coord 

390 

TYPE (neighbor IDs ) , DIMENSION  {: , :) , INTENT { IN)  ::  nID 

391 

INTEGER, INTENT (IN)  ::  ii,jj,kk 

392 

INTEGER  ::  i, j 

393 

REAL, DIMENSION  (3, 3)  ::  gridMet 

394 

REAL , DIMENS ION (3,3)  ::  jac 

395 

REAL  : :  detJac 

396 

397 

!  differencing  physical  space  (x,y,z)  W.R.T.  grid  directions  i,j,k 

398 

!  i.e.  dx/dxi  is  the  difference  of  "x"  in  the  "xi,  or  1”  coordinate 

399 

!  direction. 

400 

DO  i-1,3  !  For  each  grid  index  direction  i,j,k 

401 

DO  j-1,3  !  For  each  coordinate  direction  i.e.  x,y,z 

402 

!  Differencing  of  coordinate  'j'  in  the  * i *  direction 

403 

jac(j,i)  «  dol (coord, j,nID, i) 

404 

ENDDO 

405 

ENDDO 

406 

407 

! The  jacobianA-l  of  transformation  - 

408 

detJac  -  1 . 0d0/ (  j ac  (1 , 1 ) *  ( j ac  (2 , 2 ) * jac (3, 3) - jac (2 , 3 ) * jac  (3, 2) )  -  & 

409 

jac  (1,2) *(  jac (2, 1) * j§c(3, 3) -jac (2, 3) * jac (3, 1) )  +  & 

410 

jac(l,3)*(jac(2,l)*jac(3,2)-jac(2,2)  *  jac  (3,1)  )  ) 

411 

412 

!  These  are  the  grid  metrics - - - - 

413 

!  (1,1)  =  d_xi/d_x  ...  (2,3)  -  d_eta/d_z  ...  etc. 

414 

gridMet(l,l)  -  detJac  *  (jac (2, 2) * jac (3, 3)  -  jac (2, 3) * jac  (3, 2) ) 

415 

gridMet (1,2)  -  -detJac  *  ( jac  (1, 2 ) * jac (3, 3)  -  jac (1,  3) * jac (3, 2) ) 

416 

gridMet(l,3)  -  detJac  *  ( jac  (1, 2) * jac (2, 3)  -  jac (1, 3) ♦jac  (2, 2) ) 

417 

gridMet (2,1)  -  -detJac  *  (jac  (2, 1) * jac (3, 3)  -  jac (2, 3) * jac  (3, 1 ) ) 

418 

gridMet(2,2)  «  detJac  *  ( jac  (1, 1) * jac (3, 3)  -  jac (1, 3) * jac  (3, 1) ) 

419 

gridMet (2,3)  -  -detJac  *  ( jac  (1, 1) * jac (2, 3)  -  jac (1, 3) * jac  (2, 1) ) 

420 

gridMet (3,1)  -  detJac  *  (jac  (2, 1) *jac (3, 2)  -  jac (2, 2) * jac  (3, 1) ) 

421 

gridMet (3,2)  «  -detJac  *  ( jac  (1, 1) * jac (3, 2)  -  jac (1, 2) * jac  (3, 1 ) ) 

422 

gridMet (3,3)  =  detJac  *  (jac  (1,1) *jac(2,2)  -  jac (1, 2) * jac  (2, 1) ) 

423 

424 

END  FUNCTION  setMetrics 

425 

426 

FUNCTION  dol (vfunc , nVar , nbID, dir)  RESULT(der) 

427 

!  dol  computes  the  difference  related  to  a  first  derivative 

428 

1  gives  central  difference  for  interior  nodes  .  .  . 

429 

!  gives  forward  or  backward  difference  for  boundary  nodes  . . . 

430 

!  difference  is  used. 

431 

IMPLICIT  NONE 

432 

TYPE (func) , DIMENSION ( : ) , INTENT (IN)  ::  vfunc 

433 

TYPE (neighborlDs) , DIMENSION  (:, : ) , INTENT (IN)  : :  nbID 

434 

INTEGER, INTENT (IN)  ::  dir,nVar 

435 

REAL  ::  delta, der 

436 

437 

!  If  either  neighbor  in  the  differencing  direction  is  a  boundary  . . . 

438 

IF  (  nbID (dir, 1) %isBoundary  .OR.  nbID (dir , 2) %isBoundary  )  THEN 

439 

delta  -1.0  !  set  delta  to  1.0  (forward  or  backward  difference) 

440 

ELSE 

441 

delta  -  2.0  !  set  delta  to  2 . 0  (central  difference) 
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442 

ENDIF 

443 

!  The  computational  space  has  a  delta  of  "1.0"  grid  unit  per  index. 

444 

!  This  delta  assumption  cancels  out  when  performing  the  differencing 

445 

!  on  the  physical  variables  so  it  can  be  anything. 

446 

447 

!  This  gives  the  (ii,jj,kk)  ind 

448 

der  -  (  vfunc (nbID (dir, 2 ) %zID)  %f  (nbID (dir, 2)  %ii,  & 

449 

nbID (dir , 2 )  %  j  j  , nbID (dir, 2) %kk, nVar)  -  & 

450 

vfunc (nbID(dir,l)%zID)%f  (nbID (dir, 1)  %ii,  & 

451 

nbID (dir , 1 ) %j j , nbID (dir, 1) %kk, nVar)  )  /  delta 

452 

453 

END  FUNCTION 

454 

455 

FUNCTION  kd ( i , j )  RESULT (delta )  !  Kronnecker  delta 

456 

INTEGER,  INTENT (IN)  ::  i,j 

457 

REAL  : :  delta 

458 

459 

IF  (i  ••  j)  THEN 

460 

delta  *  l.OdO 

461 

ELSE 

4  62 

delta  *  0 . OdO 

463 

END  IF 

464 

465 

END  FUNCTION  kd 

4  66 

467 

FUNCTION  pt(i,j,k)  RESULT (epsilon)  !  permutation  tensor 

468 

INTEGER,  INTENT (IN)  ::  i,j,k 

469 

REAL  : :  epsilon 

470 

471 

IF  ( (i  —  1  .AND.  j  —  2  .AND.  k  —  3)  .OR.  & 

472 

(i  --  2  .AND.  j  —  3  -AND.  k  —  1)  .OR.  & 

473 

(i  —  3  .AND.  j  —  1  .AND.  k  —  2)  )  THEN 

474 

epsilon  *  l.OdO 

475 

ELSE  IF  ((i  -*  1  .AND.  j  »*  3  .AND.  k  —  2)  .OR.  & 

476 

(i  «  2  .AND.  j  —  1  .AND.  k  —  3)  .OR.  & 

477 

(i  --  3  .AND.  j  —  2  .AND.  k  —  1) )  THEN 

478 

epsilon  *  -l.OdO 

479 

ELSE 

480 

epsilon  *  O.OdO 

481 

END  IF 

482 

483 

END  FUNCTION  pt 

484 

485 

FUNCTION  omega (vFunc, uvwID, coord, i, j , k, nID, dir) 

486 

!  omega  computes  the  vorticity  in  the  'dir'  coordinate  direction. 

487 

!  'uvwID'  is  a  1x3  array  indicating  which  variables  are  the  velocity 

488 

!  components  in  the  x,y,z  directions.  In  the  current  data  set  this 

489 

!  should  be  set  to  uvwID  *  (/2,3,4/) 

490 

IMPLICIT  NONE 

491 

TYPE(func) , DIMENSION ( :), INTENT (IN)  ::  vFunc 

492 

TYPE(func) , DIMENS  ION ( :), INTENT ( IN)  ::  coord 

493 

INTEGER, INTENT (IN)  ::  i,j,k,dir 

494 

INTEGER, DIMENSION (:), INTENT (IN)  ::  UVWID 

495 

TYPE (neighbor IDs ) , DIMENSION  ( : , :) , INTENT (IN)  :  ;  nID 

496 

INTEGER  ; :  m,  n 

4  97 

REAL  : :  omega 

498 

499 

omega  *  0 . OdO 

500 

DO  m-1, 3 ;  DO  n-1 , 3 

501 

IF  (  pt (dir , m, n)  .EQ.  0.0  )  THEN 

502 

CYCLE 

503 

ELSE 

504 

omega  *  omega  +  & 
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Abstract  Submitted 
for  the  DFD12  Meeting  of 
The  American  Physical  Society 

Scaling  model  for  nonlinear  supersonic  jet  noise  WOUTIJN 
BAARS,  CHARLES  TINNEY,  The  University  of  Texas  at  Austin  —  Numerous 
endeavors  have  been  undertaken  to  investigate  nonlinear  propagation  of  sound  from 
jet  flows  in  range-restricted  environments.  However,  only  weak  observations  of  cu¬ 
mulative  nonlinear  effects  have  been  made  using  these  laboratory-scale  setups,  all 
the  while  being  observed  under  full-scale  conditions.  The  inconsistency  is  caused  by 
the  lack  of  rigor  in  understanding  what  the  appropriate  scaling  parameters  should 
be  for  producing  measurable  cumulative  nonlinearities  in  laboratory-scale  environ¬ 
ments.  A  scaling  model  will  be  presented  that  one  could  use  to  guide  future  studies 
aimed  at  investigating  this  unique  component  of  turbulent  mixing  noise.  At  first, 
the  important  length-scales  for  cumulative  nonlinear  waveform  distortion  -  the  shock 
formation  distance  and  the  acoustic  absorption  length  -  are  written  in  terms  of  jet 
exit  parameters.  Their  ratio,  expressed  as  the  effective  Gol’dberg  number,  is  a 
measure  of  the  strength  of  nonlinear  distortion  relative  to  that  of  dissipation.  By 
computing  the  individual  length-scales  and  this  dimensionless  ratio  for  an  experi¬ 
ment  that  is  being  designed,  one  can  estimate  the  presence  of  cumulative  nonlinear 
distortion  beforehand. 


Woutijn  Baars 
The  University  of  Texas  at  Austin 

Date  submitted:  12  Aug  2012  Electronic  form  version  1.4 


Appendix  B  -  Page  B-2 

Quarterly  Progress  Report  No.  6  -  N00014-1 1-1-0752 

Approved  for  Public  Release  -  Distribution  is  Unlimited 


Active  Control  of  Noise  from  Hot  Supersonic  Jets  |  Quarterly  Rpt.  6 


Appendix  C:  Panickar  et  al.  “Localization  of  Acoustic  Sources  in  Shock- Containing  Jet 
Flows  Using  Phased  Array  Measurements.”  AIAA  Paper  2013-0613. 


Appendix  C  -  Page  C-l 

Quarterly  Progress  Report  No.  6  -  N00014-1 1-1-0752 

Approved  for  Public  Release  -  Distribution  is  Unlimited 


Active  Control  of  Noise  from  Hot  Supersonic  Jets  |  Quarterly  Rpt.  6 


51it  A1AA  Atroipact  Sclencti  Mettlng  Including  the  New  Horizons  Forum  end  Aerospace  Exposition  Al  AA  201 3* *061 3 

07  •  10  January  2013,  Grspavlna  (Dsllas/Ft  Worth  Region),  Taxas 


Localization  of  Acoustic  Sources  in  Shock-Containing 
Jet  Flows  Using  Phased  Array  Measurements 


8 

P 

j 


1 

ll 

s 

f 

I 

z 

£ 

i 


Praveen  Panickar,*  James  P.  Erwin*  and  Neeraj  Sinha* 

Combustion  Research  and  Flow  Technology,  Inc.  (CRAFT  Tech),  Pipcrsville,  PA-18947,  USA 

Nathan  E.  Murray^  and  Gregory  Lyons^ 

National  Center  for  Physical  Acoustics ,  The  University  of  Mississippi,  University,  MS-38677,  USA 


In  this  paper,  we  shaii  examine  the  feasibility  of  using  time-resoived  hybrid  RANS-LES 
(HRLES)  simuiation  data  to  perform  noise  source  localization  studies  on  a  hot  jet  from  a 
conic  nozzle  operating  at  pressure  matched  conditions.  The  source  localization  wiii  be  per¬ 
formed  using  a  traditional  deiay-and-sum  beamforming  technique.  This  technique  requires 
time-resoived  data  on  a  phased  array  of  microphones  located  in  the  acoustic  farfieid  of  the 
Jet;  this  data  wiii  be  obtained  by  coupiing  the  HRLES  simuiation  with  a  Ffowcs  Wiiiiams 
and  Hawkings  equation  noise  prediction  code.  Using  insights  gained  from  experimentai 
beamforming,  we  shaii  show  that  beamforming  using  CFD  data  is  a  feasible,  and  poten¬ 
tially  iess  expensive  and  time-saving,  aiternative  to  constructing  complicated  phased  array 
systems  for  performing  these  calculations  on  experimentai  data. 


I.  Background  and  Motivation 

Noise  from  supersonic  jets  falls  primarily  into  two  categories,  either  turbulent  mixing  noise  or  shock  noise. 
Of  these  two,  turbulent  mixing  noise  is  comprised  of  two  sources,  both  of  which  are  broadband  in  nature: 
(a)  the  eddy  Mach  wave  radiation  which  propagates  downstream  relative  to  the  jet  flow  direction,  and  (b) 
fine  scale  turbulence  which  is  omni-directional.  Ou  the  other  hand,  shock  noise  can  be  either  narrowband 
(also  called  jet  screech)  or  broadband  (called  broadband  shock  associated  noise,  or  BSAN),  and  propagates 
towards  the  sideline  or  upstream  direction  relative,  to  the.  jet  flow.  Shock  noise  most  commonly  occurs  when 
a  convergent-divergent  nozzle  is  operated  at  off-design  conditions.  However,  jets  exhausting  from  military 
nozzles  typically  contain  shocks  even  when  operated  at  design  conditions;  this  is  primarily  because  the 
internal  contour  of  these  nozzles  are  not  shaped  in  order  to  ensure  a  smooth  expansion  of  the  gas. 

.Jet  noise  contributes  significantly  to  noise-induced  hearing  loss,  structural  degradation  of  airframes,  and 
restrictions  to  maintenance,  testing,  and  training  schedules  due  to  noise  pollution  of  communities  surrounding 
military  installations.  To  this  end,  it  is  imperative  to  gain  a  better  understanding  of  jet  noise  generation 
mechanisms  in  a  turbulent  flow,  buck  an  understanding  is  essential  if  one  is  to  construct  predictive  models 
for  jet  noise.  One  approach  that  has  been  considered  is  the  beamforming  technique  using  a  phased  array  of 
microphones.1  This  technique  has  been  shown  to  provide  an  estimate  of  the  distribution  of  noise  sources  in  a 
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flow.  Investigations  on  the  use  of  phased  arrays  in  aeroacoustics  in  general,  and  jet  noise  in  particular,  were 
conducted  by  many  researchers  including  Brooks  and  his  coworkers2^  at  NASA  Langley,  Bridges  and  his 
coworkers5'6  at  NASA  Glenn,  Papamoschou,  Morris  and  McLaughlin7  to  name  a  few.  Recently  Dougherty8 
developed  a  new  generalized  inverse  beamforming  method  for  jet  noise  and  demonstrated  its  effectiveness  in 
mapping  coherent  sources  and  extrapolating  directivity  patterns  using  a  2D  phased  array  system. 

All  the  reports  referenced  in  the  previous  paragraph  are  based  on  experimental  studies.  Working  with 
simulated  data  removes  many  of  the  encumbrances  that  are  present  while  working  in  an  experimental  facility. 
For  instance,  in  an  experimental  facility,  one  is  restricted  by  the  size  and  placement  of  microphones  that 
comprise  the  phased  array  system.  It  is  not  possible  to  use  a  phased  array  system  to  interrogate  a  cross 
plane  of  the  jet  unless  one  resorts  to  a  cage  array,  in  which  case  care  has  to.be  taken  that  the  microphones 
are  not  damaged  by  the  flow;  this  can  make  the  aperture  of  the  phased  array  larger  than  what  would  be 
practically  required  depending  on  the  frequency  of  interest.  Additionally,  moving  the  system  and  setting 
it  up  at  a  different  location  takes  a  long  time  since  each  move  requires  a  reealibration.  These  difficulties 
are  not  present  when  simulated  data  is  used  for  phased  array  calculations;  the  phased  array  can  be  located 
anywhere  in  the  flowfield.  Other  practical  limitations  present  in  experimental  testing,  such  as  the  number 
of  microphones  available  for  use,  are  also  not  present  when  looking  at  simulated  data.  The  major  limitation 
when  using  simulated  data  is  the  cost  involved  with  generating  time  resolved  data;  however,  this  difficulty 
is  gradually  being  alleviated  as  supercomputing  clusters,  with  ever-increasing  processor  counts,  become 
available,  and  by  using  these  clusters  in  parallel  to  provide  time  resolved  data  at  a  reasonably  fast  rate. 

Given  all  of  the  limitations  discussed  previously  with  experimental  phased  arrays,  it  should  be  noted  here 
that  the  phased  array  calculations  using  simulated  data  is  not  meant  to  replace,  or  eliminate,  the  need  for 
these  experiments.  Rather,  the  intent  is  to  use  the  simulated  data  to  carefully  guide  the  experiments  and 
suggest  optimized  array  configurations  that  lead  to  more  accurate  source  localization.  The.  simulated  phased 
array  data  can  be  used  in  two  specific  ways.  First,  by  using  a  large  amount  of  numerical  microphones  that 
can  be  located  anywhere  in  the  acoustic  field,  beamforming  calculations  can  be  used  in  order  to  determine 
source  characteristics  (frequency  and  location)  that  would  be  difficult  to  predict  given  the  limited  number 
of  microphones  available  to  the  experimentalist.  Secondly,  once  these  characteristics  are  known,  a  virtually 
limitless  amount  of  array  configurations  can  be  tested  using  only  the  number  of  microphones  available  for  a 
given  experiment  in  a  test  facility.  The  configurations  that  best  match  the  ideal  solution  can  be  the  focus 
of  the  experiments  for  further  testing.  The  goal  of  this  paper  is  to  present  this  methodology  and  consider 
its  application  to  shock-containing  jets  for  initial  validation,  before  expanding  out  to  other  aeroacoustics 
configurations  of  interest. 


II.  Configurations  and  Objectives 

The  current  work  involves  a  combined  experimental  and  computational  initiative  in  order  to  examine 
and  understand  the  inner  nozzle  and  near  nozzle  flow  conditions  and  their  impact  on  jet  noise.  To  this  end, 
we  shall  examine  a  nozzle  that  consists  of  a  straight  walled  convergent-divergent  section,  which  we  shall 
refer  to  as  the  conic  nozzle;  the  interested  reader  is  referred  to  Murray  et  al.9  for  nozzle  geometry  details. 
Based  on  the  throat-to-exit  area  ratio,  the  conic  nozzle  has  a  prescribed  “design”  operating  Mach  number 
of  Mj  =  1.74. 

The  primary  objective  of  the  current  work  will  be  to  investigate  the  feasibility  of  performing  phased  array 
calculations  on  simulated  jet  data.  This  feasibility  will  be  gauged  by  performing  beamforming  calculations  on 
the  conic  nozzle  with  the  centerbody  operating  at  its  design  (or  pressure- matched)  condition.  This  primary 
objective  can  be  broken  up  into  the  following  sub-objectives: 
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1.  Compare  overall  and  spectral  acoustic  characteristics  of  experimental  and  computational  configurations 
in  the  farfield  in  order  to  ensure  agreement  between  the.  two. 

2.  Acquire  fluctuating  pressure  data  using  a  linear  phased  array,  for  the  conic  nozzle  configuration  with 
the  ccnterbody  operating  at  pressure  matched  conditions,  in  the  anechoic  jet  laboratory  at  the  NCPA 
in  order  to  determine  an  acceptable  number  of  ensembles  for  constructing  the  cross-spectral  matrix. 

3.  Perform  high-fidelity  HRLES  simulations  of  this  configuration,  and  obtain  adequate  time-resolved  data 
on  the  ADS  (as  dictated  by  the  results  of  the  previous  item). 

4.  Using  FW-H  equations,  propagate  the  acoustics  on  the  ADS  to  microphone  locations  on  the  linear 
array  that  was  used  in  the  experimental  measurements. 

5.  Perform  beamforming  calculations  of  the  experimental  and  computational  configurations  and  compare 
the  source  localization  results  provided  by  both  methods. 

The  experiments  and  computational  measurements  will  be  performed  on  heated  jets  with  stagnation  tem¬ 
perature  of  1005  K  (1350  °F)  unless  specified  otherwise. 

III.  Computational  Capabilities  and  Experimental  Setup 

The  computational  simulations  were  performed  using  CRAFT  Tech’s  CRAFT  CFD®  LES  code10, 1 1  which 
is  a  structured  Navicr-Stokes  solver  that  uses  an  upwind  differencing  scheme  that  is  fifth  order  accurate  in 
space  and  various  time  integration  schemes  including  a  second  order  implicit  alternating  direction  implicit 
(ADI)  scheme  and  various  explicit  Runge-Kutta  schemes.  The  code  is  highly  parallelized  using  domain 
decomposition,  with  linear  speed-up  parallel  performance  on  several  large  supercomputers  for  very  large 
grids.12  The  calculations  were  performed  using  a  hybrid  RANS-LES  (HRLES)  model.13  The  primary 
advantage  of  using  HRLES  is  the  ability  to  include  the  internal  nozzle  in  the  jet  noise  simulations.  By 
including  the  internal  nozzle  it  becomes  possible  to  smoothly  transition  the  subgrid  eddy  viscosity  from 
within  the  nozzle  boundary  layer  to  the  LES.  The  intent  of  the  HRLES  simulations  is  to  more  accurately 
describe  the  nozzle  and  LES  interface  boundary  condition  instead  of  dialing  in  a  steady  boundary  condition 
from  a  precursor  RANS  solution  (zonal  RANS/LES).  The  HRLES  framework  allows  for  feedback  between 
the  RANS  and  LES  regions  and  can  help  appropriately  capture  realistic  fluctuations  inside  of  the  nozzle 
given  sufficient,  and  often  times  restrictive,  grid  resolution.  Currently  the  eddy  viscosity  predicted  by  the 
HRLES  model  tends  to  smooth  the  initial  shear  layer  and  negatively  impact  the  noise  levels  in  the  sideline 
and  upstream  directions.  In  the  downstream  direction,  the  primary  component  of  noise  is  due  to  turbulent 
mixing  and  this  component  is  captured  with  sufficient  accuracy  by  the  zonal  RANS/LES  method.  Thus,  by 
using  HRLES,  one  can  expect  to  improve  noise  predictions  in  the  sideline  and  upstream  directions  and  show 
better  agreement  with  experiments.  The  computational  grid  for  the  current  initiative  contains  approximately 
7  million  cells  and  extends  60 Rj  (Rj  is  the  nozzle  radius)  in  the  axial  (X)  direction  and  had  a  radius  of 
30Rj  at  its  widest  extent.  Downstream  of  this,  there  is  a  buffer  zone  that  extends  a  further  30 Rj  in  order  to 
damp  out  unwanted  reflections  at  the  LES  flow  boundaries.  Inside  of  the  buffer  zone,  rapid  grid  stretching 
is  employed  along  with  a  first  order  upwind  scheme. 

Many  researchers  have  contributed  to  the  development  of  the  CFD  acoustic  analogy  method  for  jet 
noise.14-17  In  this  paper  we  shall  make  use  of  the  Ffowcs- Williams  and  Hawkings  (FW-H)  equation  to 
compute  the  far-field  noise  at  desired  microphone  locations.  This  solution  consists  of  an  integration  over  a 
surface  surrounding  the  noise  generating  mechanisms  of  the  flow. 18,19  In  order  to  perform  this  integration, 
flow  variables  are  recorded  on  this  “permeable”  acoustic  data  surface  (ADS)  with  the  desired  time  resolution. 
Figure  1  shows  the  HRLES  computational  domain  with  the  ADS  in  place.  The  time  step  for  the  CFD 
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Figur*  1.  Computation*!  domain  with  acoustic  data  surfaca  (ADS),  shown  shaded,  used  to  predict  farfl  Jd  nolsa  using 
tha  Ffowcs-Wlilimms  and  Hawkings  (FW-H)  equations. 


simulations  in  the  current  effort  was  selected  to  be  5.0e-8  secouds.  For  farfield  uoise  and  beanifonniug 
calculations  a  sample  rate  of  200  kHz  was  determined  to  be  satisfactory,  and  acoustic  variables  were  recorded 
on  the  ADS  every  100  LES  tiraesteps.  The  LES/FW-H  codes  described  here  have  been  extensively  used  for 
the  prediction  of  jet  noise  for  a  wide  variety  of  subsonic  and  supersonic  jets  in  both  free  and  impinging  jet 
configurations.20-  22 

In  order  to  supplement  the  CFD  computations,  some  experimental  phased  array  tests  were  conducted  in 
the  Anechoic  Jet  Laboratory  (AJL)  at  the  University  of  Mississippis  National  Center  for  Physical  Acoustics; 
facility  details  can  be  found  in  Murray  et  al.9  This  phased  array  system  consists  of  32  exponentially  spaced 
Kulitc  pressure  transducers  (type  XT-140,  dynamic  range  m  0-100  psia),  mounted  on  a  stainless  steel  rod. 
This  array  was  positioned  G0.9G  cm  (24  inches)  below  the  nozzle  exit  and  offset  7.G4  cm  (2.94  inches)  in  the 
axial  direction;  the  downstream  end  of  the  rod  was  inclined  at  an  angle  of  15°  to  the  nozzle  axis.  Figure 
2  shows  an  image  of  the  experimental  setup  with  the  nozzle  and  the  pressure  transducers  mounted  in  the 
rod.  Digitized  fluctuating  pressure  data  from  the  array  was  recorded  and  stored  on  a  National  Instruments 
PXIe  (PCI  extensions  for  Instrumentation-express)  system.  The  PXIe  system  has  four  PXIe-4331  cards 
(8  channels,  24  bits  of  resolution,  102.4  kHz  maximum  sample  rate).  Each  acquisition  recorded  1048576 
samples  per  channel  at  100  kHz.  For  a  given  configuration,  at  least  10  such  acquisitions  of  data  were 
recorded.  Spectral  data  were  calculated  by  dividing  the  data  in  a  given  run  into  blocks  of  8192  points  each; 
this  yielded  a  frequency  resolution  of  12.21  Hz.  Depeudiug  on  the  number  of  files  used  to  do  the  spectral 
calculations,  the  number  of  ensembles  that  could  be  averaged  for  the  cross  spectral  matrix  ranged  from  1 
(only  part  of  one  file  used)  to  1280  (all  ten  files  used). 

IV.  Results  and  Discussion 

Figure  3  shows  the  LES  time-averaged  stream  wise  velocity  component,  £/,  of  the  conic  nozzle  operating 
at  pressure  matched  conditions.  The  sharp  throat  and  the  conic  diverging  section  (as  opposed  to  a  smooth 
throat  and  divergent  section  that  would  be  produced  using  the  method  of  characteristics)  results  in  two 
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Figure  2.  Experimental  phased  array  setup. 


separate  shock-trains;  the  first  originating  at  the  sharp  throat,  and  the  second  originating  at  the  nozzle  exit. 
Both  shock-trains  propagate  downstream  along  the  jet  column  as  seen  in  the  figure.  In  addition  to  the  shock- 
train,  one  can  also  see  the  slip  line  along  the  jet  axis  where  the  oblique  shock-cells  coalesce.  Figure  4  shows 


Figure  S.  Mean  straamwlaa  velocity  contour  of  the  conic  noxxle  with  upstreem  centerbody  operetlng  et  pressure 
matched  conditions. 


the  normalized  (using  the  jet  exit  velocity)  streamwise  velocity  along  the  jet  axis.  Tills  figure  clearly  shows 
the  oscillations  in  the  mean  velocity,  which  are  a  result  of  the  jet  undergoing  rapid  successive  expansions 
and  contractions,  which  shows  up  in  the  flow  visualization  as  shock-cell  structures. 
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Figure  4.  Normalized  mean  atreamwlze  velocity  along  tha  centerline  of  the  jet. 


IV. A.  Comparison  of  Farfield  Acoustic  Characteristics 

Prior  to  performing  phased  array  calculations  on  the  simulated  jet  data,  it  is  important  to  compare  the 
simulation  and  experiments  in  order  to  ensure  agreement  between  the  two.  In  the  current  paper,  this 
agreement  will  be  verified  by  comparing  the  farfield  acoustic  data  measured  using  the  arc  array  of  microphones 
in  the  AJL  to  the  acoustic  levels  at  the  same  microphone  locations  predicted  using  LES/FW-H  methodology. 
Figure  5  shows  a  comparison  of  the  overall  sound  pressure  levels  (OASPL)  between  the  experimental  data 
and  LES/FW-H  calculations  at  the  various  azimuthal  location  which  are  measured  with  reference  to  the 
upstream  location.  It  can  be  seen  that  the  experiments  and  simulations  agree  to  within  ±2  dB  over  most 
of  the  measurement  angles.  The  peak  level  and  directivity  agrees  well  with  the  experimental  measurements, 
including  the  OASPL  directivity  at  upstream  angles,  which  is  a  characteristic  of  shock-containing  jets.  The 
agreement  worsens  at  the  downstream  angles,  but  this  can  be  attributed  to  the  coarsening  of  the  HRLES 
grid  at  the  downstream  locations  and  the  finite  length  of  the  ADS  (extending  only  to  60 Rj). 

Figure  6  shows  a  comparison  of  the  noise  spectra  at  three  different  farfield  microphone  locations:  135° 
(turbulent  mixing  noise  direction),  90°  (sideline  direction),  and  60°  (shock  noise  direction).  Once  again,  the 
agreement  between  the  experimental  measurements  and  simulations  are  acceptable  to  within  the  restrictions 
imposed  by  the  grid  resolution  (~  10  kHz).  In  the  downstream  direction,  the  Mach  wave  radiation  peak 
agrees  very  well  with  the  experimental  measurements.  The  over-prediction  of  low-frequency  noise  at  all 
angles  is  due  to  the  finite  LES  runtime  before  noise  sampling,  and  the  existence  of  a  long-scale  acoustic 
instability  that  exists  between  the  inflow  and  outflow  of  the  LES  domain.  In  the  sideline  and  upstream 
directions,  there  appears  to  be  a  build-up  of  energy  just  before  the  grid  cut-off  frequency  of  around  10 
kHz.  This  build-up  is  potentially  due  to  the  over-prediction  of  eddy-viscosity  by  the  HRLES  model,  and 
holding  the  initial  shear  layer  steady  for  too  far  downstream.  When  the  shear  layer  destabilizes,  the  grid 
resolution  in  this  region  is  fine  and  therefore  leads  to  the  over-prediction  of  mid-to-high  frequency  noise  just 
before  the  cut-off.  To  improve  the  sideline  and  upstream  noise  predictions,  more  grid  resolution  is  required. 
Nevertheless,  the  coarse  grid  used  in  this  study  yields  acceptable  noise  predictions  which  will  help  guide  the 
development  of  the  source  localization  in  the  remainder  of  the  paper. 
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Flgura  B.  Comparison  of  tha  overall  sound  pressure  levels  between  experimental  measurements  and  FW-H  calculations 
of  the  simulated  Jet. 


IV. B.  Experimental  Phased  Array  Measurements 

Prior  to  performing  phased  array  calculations  on  the  HRLES  data,  a  rudimentary  set  of  experiments  were 
conducted  on  the  conic  nozzle  configuration  with  the  centerbody  using  the  experimental  linear  phased  array 
system.  The  nozzles  were  operated  at  cold  and  heated  conditions  and  fluctuating  pressure  data  was  recorded 
by  the  linear  phased  array  system.  The  primary  objective  of  these  experiments  was  to  provide  guidance  to 
the  CFD  simulations  regarding  an  acceptable  number  of  ensembles  required  to  generate  the  cross-spectral 
matrices  used  for  beamforming  calculations.  Tb  this  end,  fluctuating  pressure  data  was  acquired  by  the 
linear  array  for  the  cold  nozzle  (To  =  Tom^ent)  operated  at  overexpanded  conditions  (Mj  =  1.55).  Figures  7 
show  sample  spectra  and  cross-correlations  coefficients  calculated  for  the  first  six  microphones  in  the  array 
for  this  cold  jet  operated  at  overexpanded  conditions.  The  spectra  and  cross-correlations  shown  in  figure 
7  were  calculated  using  128  ensemble  averages.  The  spectra  clearly  show  the  presence  of  both  shock  noise 
components,  i.e.,  screech  and  broadband  shock  associated  noise.  The  screech  tone  can  be  clearly  seen  at 
approximately  2.6  kHz  and  the  broadband  shock  noise  hump  is  seen  at  higher  frequencies  between  4  and  5 
kHz. 

In  contrast,  as  shown  in  figure  8,  for  heated  jets  operating  at  the  same  overexpanded  pressure,  screech 
tones  can  no  longer  be  clearly  discerned.  Additionally,  it  can  be  seen  that  the  broadband  hump  has  moved  to 
higher  frequencies  and  the  correlation  levels  between  the  sensors  has  reduced.  This  points  to  a  fundamental 
change  in  the  acoustic  sources  that  we  shall  examine  using  phased  array  calculations. 

The  data  from  the  phased  array  system  is  reduced  using  the  classical  delay-and-sum  (DAS)  beamforming 
method  in  the  frequency  domain.  This  method  assumes  that  the  noise  sources  are  incoherent,  point  sources. 
In  order  to  arrive  at  the  bcamform  map,  the  following  procedure  is  followed: 

*  Construct  the  cross-correlation  matrix  for  each  narrowband  frequency,  S /  (n,m),  where  n  and  m  are 
the  microphone  indices,  and  /  is  the  narrowband  frequency;  for  the  32  microphones  used  in  this  study, 
this  matrix  has  size  32  x  32  for  each  frequency. 


7  of  17 


American  Institute  of  Aeronautics  and  Astronautics 


Appendix  C  -  Page  C-8 

Quarterly  Progress  Report  No.  6  -  N00014-11-1-0752 

Approved  for  Public  Release  -  Distribution  is  Unlimited 


Active  Control  of  Noise  from  Hot  Supersonic  Jets  [  Quarterly  Rpt.  6 


I 

a 


r 

i 

>. 

i 

§ 


! 

u. 

g 

i 

z 


Conic  nozzle  with  centertxxfy  -  pressure  matched  operation 


Flgura  6.  Comparison  of  th«  farflald  spactral  charactaristica  batwaan  sxpsrimcntsl  mMiurctncnti  and  HRLES  ilmuU* 
tlons  for  ths  conic  nomsU  with  cantarbody  oparatlng  at  praasura  matchad  conditions. 


8  of  17 


American  Institute  of  Aeronautics  and  Astronautics 


Appendix  C  -  Page  C-9 

Quarterly  Progress  Report  No.  6  -  N00014-11-1-0752 

Approved  for  Public  Release  -  Distribution  is  Unlimited 


Active  Control  of  Noise  from  Hot  Supersonic  Jets  [  Quarterly  Rpt.  6 


Frequency  [Hz] 


i-T 

i 

Uu 

s 

f 


(•)  Spectra 


(b)  Cross  correlation  coefficients 

Figure  7.  (•)  Spectre  end  (b)  cross  correietion  coefficients  for  the  first  six  sensors  In  the  experiments!  lineer  phased 

array.  Cold  Jet,  overexpanded  operation,  Mj  *=  1.55. 


•  Delete  the  diagonal  element  of  the  cross-correlation  matrix  at  each  frequency  to  give  the.  modified 
cross-spectral  matrix  S/  (m,  n). 


•  The  bcamforming  result  is  calculated  using  the  expression: 


ilk  =  Wk§Wk, 


(1) 


where  u>*  represents  the  steering  vector  for  the  microphones  to  a  steering  point  A:,  and  is  expressed  as: 


w 


e  c 


\Xm  -  Xk\  ’ 


(2) 


where  m  is  the  microphone  index,  w  =  2nf  is  the  circular  frequency,  c  is  the  speed  of  sound,  and  x 
is  the  vector  coordinate  location,  and  |£m  -  £*|  represents  the  distance  between  the  mth  microphone 
location  and  the  kth  scanning  plane  location. 
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(b)  Cross  correlation  coefficients 

Figure  8.  (•)  Spectra  and  (b)  cross  correlation  coefficients  for  tha  first  six  sensors  In  the  experimental  linear  phased 

array.  Heated  Jat,  Tj  «  1005  K,  pressure-matched  operation,  m  1.74. 


Using  this  method  the  beamforming  result  is  calculated  on  a  plane  of  interest,  in  this  case  the  symmetry 
plane  (Z  =  0),  oriented  in  the  stream  wise- vertical  (X,y)  direction.  Following  this,  the  beamforming  results 
are  presented  as  contour  plots  corresponding  to  a  specific  1  /3rd  octave  band  of  interest.  This  is  done  by 
summing  the  beamforming  results  at  each  narrowband  frequency  comprising  the  1  /3rd  octave  band. 

For  the  linear  phased  array  system,  the  beamforming  output  is  on  a  scanning  line,  that  extends  along 
the  nozzle  axis,  rather  than  on  a  scanning  plane.  Results  are  then  presented  as  contour  plots  on  the  axial 
location- frequency,  ( X}f)1  plane.  For  phased  array  calculations  using  CFD  data,  results  are  presented  in 
the  traditional  manner,  on  the  ( X ,  T)  plane. 


IV.B.l.  Determining  Optimal  Number  of  Ensembles 

When  dealing  with  time-resolved  HRLES  simulation  data,  it  must  be  remembered  that  it  is  generally  not 
feasible  to  acquire  as  much  time-resolved  data  as  one  can  in  an  experimental  study.  For  instance,  the 
time-resolved  experimental  data  was  acquired  over  a  period  of  ~  10  seconds;  to  generate  the  same  amount 
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of  simulation  data  would  take  a  prohibitive  amount  of  computational  time.  For  this  reason,  a  study  was 
undertaken  in  order  to  determine  the  amount  of  time-resolved  data  required  for  acceptable  beamforming 
results.  This  study  was  performed  on  the  cold  conic  nozzle  with  the  centerbody  operating  at  overexpanded 
conditions  (spectra  and  cross-correlations  shown  in  figure  7)  since  this  configuration  had  distinct  frequency 
regimes  of  interest  (screech  tones  and  broadband  shock-associated  noise)  that  could  be  examined  for  the 
purpose  of  this  study. 

Figures  9(a)  and  9(b)  show  the  delay-and-sum  beamforming  results  using  128  ensembles  for  the  cold, 
overexpanded  conic  nozzle  configuration.  This  figure  clearly  shows  the  broadband  shock  associated  noise, 
between  4-6  kHz,  is  localized  within  10  inches  of  the  nozzle  exit.  Additionally,  the  screech  tones  at  2.6 
kHz,  seen  better  in  figure  9(b),  extend  intermittently  (corresponding  to  the  shock-cell  intcrmittcncy)  for 
almost  the  entire  axial  length  scanned.  Reducing  the  number  of  ensembles  to  8  reduces  the  resolution  of  the 
beamforming  results  significantly  as  shown  in  figures  9(c)  and  9(d);  however,  one  is  still  able  to  distinguish 
the  broadband  shock- associated  uoise  and  screech  tone  components  from  the  noise  fioor.  Further  reduction  of 
the  number  of  ensembles  to  4,  as  shown  in  figures  9(c)  and  9(f),  deteriorates  the  quality  of  the  bcamforming 
results  to  where  it  becomes  unacceptable.  For  this  reason,  it  was  decided  to  run  the  HRLES  simulations  for 
enough  iterations  that  would  provide  6  ensembles  for  averaging.  In  order  to  run  the  time-resolved  iterations 
of  the  HRLES  simulations  we  have  used  the  High  Performance  Computing  (HPC)  resources  provided  by 
the  Department  of  Defence  (DoD).  It  must  be  noted  that  it  takes  approximately  250,000  CPU  hours  to  run 
the  required  number  of  iterations,  which  corresponds  to  approximately  55  days  of  continuous  run  time  when 
using  190  CPUs  operating  in  parallel. 

IV.C.  Phased  Array  Calculations  Using  HRLES  Data 

For  the  CFD  data,  the  phased  array  calculations  were  conducted  using  two  different  array  designs.  Both 
designs  used  a  2D  array  of  microphones  that  were  located  parallel  to  the  streamwise- vertical,  (X,K),  plane 
of  the  jet  and  offset  by  a  distance  of  1.83  m  (72  inches)  in  the  Z  direction.  The  first  design  was  similar  to  the 
JEDA  design  discussed  in  Brooks  et  al.4  and  consisted  of  41  microphones  arranged  in  5  concentric  circles 
each  containing  8  microphones  and  an  additional  microphone  at  the  center.  The  second  design  consisted 
of  45  microphones  arranged  in  5  logarithmic  spiral  arms,  each  containing  9  microphones.  Figure  10  shows 
the  microphone  locations  of  both  phased  array  measurement  systems.  As  seen  from  the  figure,  both  phased 
array  systems  have  the  same  effective  aperture  and  size. 

As  mentioned  previously,  traditional  delay-and-sum  beamforming  calculations  with  diagonal  removal  of 
the  cross  spectral  matrix  were  performed,  using  both  phased  array  systems  at  various  l/3rd  octave  band 
center  frequencies.  The  scanning  plane  for  these  calculations  was  the  (Xy  Y)  plane  at  Z  —  0.  For  reference, 
the  nozzle  exit  plane  corresponds  to  X  =  0.  The  results  of  these  calculations  are  presented  as  stacks  of 
bcamform  maps:  figure  11  shows  the  bcamform  maps  for  the  phased  array  system  consisting  of  microphones 
arranged  in  concentric  circles  and  figure  12  shows  the  beamform  maps  for  the  phased  array  system  with  the 
log-spiral  arrangement  of  microphones.  These  figures  show  that  the  dominant  noise  source  for  the  jet  is 
located  downstream  of  the  nozzle  exit  at  the  lower  l/3rd  octave  band  center  frequencies.  For  higher  center 
frequencies,  the  source  regions  are  located  closer  to  the  nozzle  exit  and  are  smaller  in  size  compared  to  those 
at  the  lower  center  frequencies.  This  result  is  in  broad  agreement  with  experimental  work  reported  in  Brooks 
et  al.4  and  Podboy  et  al.6  Moreover,  it  is  seen  that  the  phased  array  system  using  the  logarithmic  spiral 
arrangement  of  microphones  does  not  show  striations  in  the  beamforming  results  similar  to  what  is  seen 
for  the  phased  array  system  using  the  concentric  circular  array  of  microphones;  this  is  especially  evident  at 
higher  center  frequencies. 

Finally,  for  the  purposes  of  comparison,  we  present  beamforming  results  from  a  third  array  design.  This 
array  is  similar  to  the  logarithmic  spiral  array  shown  in  figure  10  and  uses  the  same  number  of  microphones. 
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Figure  9.  Beamforming  results  using  the  linear  array  of  transducers  on  a  cold  conic  nossle  with  the  centerbody  operating 
at  pressure  matched  conditions. 
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Figure  10.  Microphone  locations  for  the  phased  array  systems  used  In  the  current  study.  The  microphones  on  tha 
concentric  phased  array  are  shown  as  filled  black  circles,  and  the  microphones  on  the  log-splrel  array  ere  shown  as  filled 
blue  circles. 


However,  this  array  is  ~  1.5  times  larger  in  diameter;  henceforth,  we  shall  refer  to  the  spiral  array  shown  in 
figure  10  as  the  small  aperture  spiral  array  (SASA)  and  this  new  design  as  the  large  aperture  spiral  array 
(LAS  A).  The  beamforming  results  using  the  LAS  A  are  shown  in  figure  13.  This  figure  shows  the  merit  of 
using  a  phased  array  of  larger  aperture  for  beamforming  calculations  at  the  lower  eenter  frequencies.  It  can 
be  seen  that  for  the  3.15  kHz,  5  kHz  and  8  kHz  center  frequencies,  the  beamforming  results  using  the  LASA, 
figure  13,  have  better  resolution  than  their  counterparts  in  figure  12  that  were  calculated  using  the  SASA. 
However,  the  disadvantage  of  using  the  LASA  ean  be  seen  from  the  beamforming  results  at  higher  eenter 
frequencies  where  the  resolution  is  diminished  when  compared  to  the  SASA.  This  reduction  in  resolution 
at  higher  frequencies  can  be  attributed  to  spatial  aliasing  in  the  LASA,  since  the  microphones  are  located 
relatively  further  apart  as  compared  to  the  SASA. 

Beamforming  calculations  using  simulated  data  can  thus  be  reliably  used  to  determine  acoustic  source 
locations  in  a  jet  flow.  Due  to  the  quick  turnaround  time  for  these  calculations,  a  large  number  of  phased 
array  designs  can  be  tested  in  a  relatively  short  period  of  time.  This  method  has  the  potential  for  providing 
guidance  for  optimal  phased  array  designs  and  locations,  which  can  then  be  implemented  and  verified  in  an 
experimental  setup. 


V,  Concluding  Remarks 

The  beamforming  results  presented  in  this  paper  serve  to  underscore  the  feasibility  of  using  simulated 
data  in  conjunction  with  a  noise  prediction  code  for  the  purposes  of  identifying  acoustic  sources  In  a  flow. 
It  is  hoped  that  these  studies  ean  be  used  to  provide  insight  into  acoustic  characteristics  of  various  flows 
of  interest  to  the  aeroacoustics  community  prior  to  manufacturing  and  testing  expensive  models  for  the 
laboratory.  In  this  sense,  the  beamforming  calculations  using  simulated  data  have  been  conceptualized  in 
order  to  provide  guidance  for  optimal  array  designs  for  experimental  efforts.  Simulations  of  various  other  jet 
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flow  configurations  of  interest,  as  well  as  various  other  aeroacoustic  flows,  are  currently  underway.  We  hope 
to  present  results  from  these  studies  at  future  conferences. 
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Figure  13.  B&e informing  maps  for  the  jet  from  the  conic  nommie  with  the  centerbody  operated  at  pressure  matched 
condition  using  the  lerge  eperture  logarithmic  spiral  microphone  array. 
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•  Conclusions  and  future  work. 
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Motivation  and  Problem  Definition 


•  Jet  noise  from  military  aircraft  engines  fall  into  two  categories:  turbulent  mixing 
noise  and  shock  noise, 

•  Better  understanding  of  the  characteristics  of  these  components  is  necessary  to 
further  efforts  in  development  of  control  and  suppression  strategies, 

-  Changes  with  jet  operating  conditions  (NPR  and  TTR), 

-  Changes  with  nozzle  geometry, 

•  Beamforming  techniques  can  provide  valuable  insight  into  the  location  and 
distribution  of  acoustic  sources  in  jet  flows, 

•  These  techniques,  combined  with  high  fidelity  simulation  data,  can  be  used  to 
guide  the  design  of  acoustic  suppression  devices  in  a  more  cost  effective  way 
compared  to  expensive  fabrication  and  laboratory  testing. 

Perform  beamforming  calculations  on  simulated  jet  flows  in  order  to  demonstrate 
the  feasibility  of  obtaining  acoustic  source  localization  and  distribution 
information ;  identify  important  parameters  and  establish  the  guidelines  necessary 
for  reliably  performing  these  calculations. 
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Flow  Configuration  and  Geometry 


•  2  inch  exit  diameter  CD  nozzle  with  conic  convergent  and  divergent  sections, 
operated  with  a  centerbody  upstream  of  the  convergent  section, 

•  Nozzle  design  Mach  number  based  on  area  ratio  =  1.74, 

Nozzle  stagnation  temperature  =  1 350  °F. 


•  Nozzle  stagn 
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CFD  Methodolgy 


•  ~  7  million  grid  points, 

•  CRAFT  CFD®  solver  using  hybrid  RANS-LES  (HRLES)  model  with  time  step 
=  5e-8  seconds, 

•  5th  order  (spatially)  upwind  differencing  scheme  and  various  time  integration 
schemes  including  a  2nd  order  ADI  and  various  Runge-Kutta  schemes, 

•  Parallelized  using  MPI. 

. "r .  1 
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CFD  Methodology 


HRLES  Simulation  -  Farfield  Acoustics 


•  Time  resolved  data  recorded  on  an  acoustic  data  surface;  this  data  is  used 
to  compute  the  noise  at  desired  farfield  microphone  locations  using  the 
Ffowcs- Williams  and  Hawkings  (FW-H)  method. 
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Experimental  Facility 


Experimental  Facility 
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•  Time  resolved  simulations:  CFD  timestep  =  5.0e-8  seconds, 

•  Wake  from  centerbody  clearly  seen 

•  Double-diamond  shock  cell  structure  due  to  shocks  downstream  of  the  sharp 
throat  and  at  the  nozzle  exit  can  also  be  clearly  seen. 
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HRLES  Simulation  -  Spectra 


Conic  nozzle  with  centerbody  -  pressure  matched  operation 
130( - ! — ■■  - - — ■  ■■■■■■, - - - ~| 
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HRLES  Simulation  -  Spectra 


Conic  nozzle  with  centerbody  -  pressure  matched  operation 
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HRLES  Simulation  -  Spectra 
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Phased  Array  Calculations 


Pros  and  Cons  of  Phased  Array  Calculation  in  CFD  Domain 

•  Use  acoustic  data  surface  and  Ffowcs-Williams  Hawking  method  to  calculate 
time  resolved  data  at  desired  microphone  location;  these  locations  may  not  be 
feasible  in  a  laboratory  setting  due  to  spatial  or  flowfield  restrictions. 

•  Allows  for  rapidly  reconfigurable  phased  arrays  of  different  designs  that  can  be 
used  to  simultaneously  “view”  the  jet  flowfield, 

-  Can  use  different  phased  array  design  for  different  frequencies  of  interest 
without  associated  costs  of  building  and  setting  up  a  new  array. 

•  Allows  for  use  of  many  more  microphones  than  practically  feasible  in 
experimental  work. 

•  The  challenge  lies  in  obtaining  ‘enough’  ensembles  in  order  to  sufficiently 

resolve  cross-correlation  estimates  between  the  microphones.  j 
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Experimental  Phased  Array  Study 

Effect  of  number  of  ensembles 

•  Linear  array  of  32  microphones  with  logarithmically  varying  spacing  between 
them, 

•  Selected  configuration  of  cold,  overexpanded  jet  ( Mj  =  1.55)  from  a  conic  nozzle 
with  a  centerbody, 

•  Sample  rate  =  100  kHz,  1024x  1024  samples  acquired. 

•  Traditional  delay-and-sum  (DAS)  beamforming  calculations. 
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Experimental  Phased  Array  Study 

128  Ensembles 

•  Blocksize  of  8192  points,  A /  =  12.21  Hz... 
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Experimental  Phased  Array  Study 


8  Ensembles 

•  Blocksize  of  4096  points,  A /  =  24.41  Hz... 
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Experimental  Phased  Array  Study 
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•  Prefer  between  4-8  ensembles. 
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4  Ensembles 

Blocksize  of  4096  points,  A/  **  48.83  Hz... 


Phased  Array  Design 


•  Two  array  orientations  considered:  in  shock-noise  direction  and  in  peak  turbulent 
mixing  noise  direction. 

•  Two  types  of  array  design  considered: 

-  Microphones  arranged  in  concentric  circles  (8  microphones  per  arm,  5  arms 
+  one  microphone  at  center *  *  total  41  microphones), 

-  2D  log-spiral  arrangement  of  microphones  (9  microphones  per  arm,  5  arms 
=  total  45  microphones), 
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CFD  Phased  Array  Methodology 


•  Using  FWH  method,  obtain  time  resolved,  fluctuating  pressure  data  at  phased 
array  microphone  locations, 

•  Calculate  and  save  cross-spectral  matrices  (CSM)  at  each  narrowband  frequency, 
G f  (m,  n),  where  m  and  n  are  microphone  indices  and  /  is  the  narrowband 
frequency, 

-  Hanning  window  applied  in  time  domain, 

-  6  bins  containing  4096  points  each,  frequency  resolution,  A /  =  48.83  Hz, 

-  Diagonal  elements  of  the  CSM  are  deleted  (reduces  self  noise  from 
microphones), 

•  Perform  classical  beamforming  (a.k.a  delay-and-sum  beamforming)  calculations 
in  the  X  —  Y  (streamwise-vertical)  plane  with  spatial  resolution  of  Ax  =  Ay  = 
0.25  inches. 

•  Array  offset  for  2D  phased  arrays  =  72  inches. 

•  Beamforming  results  presented  in  l/3rd  octave  bands  by  summing  results  at  each 
constituent  narrowband  frequency  between  the  high  and  low  frequency  limits  of 
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CFD  Octave  Band  Beamforming  Results 


2D  Circular  Array 
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CFD  Octave  Band  Beamforming  Results 


2D  Spiral  Array 


CFD  Octave  Band  Beamforming  Results 


2D  Circular  Array 
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CFD  Octave  Band  Beamforming  Results 


2D  Spiral  Array 
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Insights  from  Preliminary  Beamforming 
Calculations 

•  DAS  beamforming  shows  the  presence  of  acoustic  sources  that  appear  to  move 
closer  to  the  jet  exit  as  the  l/3rd  octave  band  center  frequency  is  increased. 

-  At  the  lower  frequency  the  source  region  extends  as  much  as  12-15  nozzle 
diameters  downstream  of  the  nozzle  exit, 

-  At  the  higher  frequencies,  the  sources  appear  as  distinct  packets,  probably 
corresponding  to  the  shock  cell  pattern. 

•  Results  emphasize  the  feasibility  of  beamforming  calculations  on  CFD 
simulations  as  a  viable  technique  to  study  acoustic  source  characteristics  of 
various  jet  configurations. 

•  Results  are  in  qualitative  agreement  with  references  showing  experimental 
beamforming  studies  on  hot,  shock-containing  jets  (Dougherty  and  Podboy 
(AIAA  Paper  2009-3186),  Brooks,  Humphreys  and  Plassman  (AIAA  Paper 
2010-3780),  and  Podboy,  Bridges  and  Henderson  (2010  NASA  TM)). 
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Effect  of  Array  Size 

Small  and  Large  Aperture  Spiral  Arrays 
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Effect  of  Array  Size 


Small  Aperture  Spiral  Array 
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Effect  of  Array  Size 


Large  Aperture  Spiral  Array 


Effect  of  Array  Orientation 


Large  Aperture  Spiral  Array,  Shock  Noise  Direction 
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Effect  of  Array  Orientation 


Conclusions  and  Future  Work 


•  Beamforming  of  simulated  data  is  a  feasible  alternative  to  building  experimental 
phased  arrays  for  localizing  acoustic  sources  in  a  jet  flow 

-  Eliminates  restrictions  on  number  of  sensors  that  can  be  used, 

-  Allows  for  design  and  calculations  using  multiple  sensor  arrays  that  would 
not  be  practically  feasible  (space  and  cost  considerations), 

•  Produce  simulation  data  with  higher  frequency  resolution  to  resolve  locations  of 
high-frequency  sources, 

•  Extend  beamforming  methods  to  include  different  free  jet  flow  configurations 
and  geometries, 

•  Determine  dependance,  if  any,  of  amount  of  time-resolved  data  required  given 
different  flow  situations, 

•  Incorporate  additional  beamforming  methods  (DAMAS,  GINV  and 
CLEAN-SC). 
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Application  of  MHz  frame  rate,  high  dynamic  range 
PIV  to  a  high-temperature,  shock-containing  jet 
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An  experimental  investigation  is  described  for  the  study  of  a  high-temperature,  shock- 
containing  jet  emanating  from  a  conic-section,  converging-diverging  nozzle.  Particie  image 
veiocimetry  (PIV)  measurements  were  acquired  for  a  flow  fleid  centered  axially  at  the 
end  of  the  jet  potential  core  and  radially  along  the  lower  half  of  the  shear  iayer.  For  aii 
cases  the  nozzle  was  operated  at  over-expanded  conditions,  and  PIV  images  were  acquired 
through  the  combined  use  of  a  pulse  burst  iaser  and  a  high-speed,  gated  intensified  CCD 
framing  camera.  Because  the  system  could  acquire  sequences  of  16  images  at  MHz  frame 
rates,  temporally  resolved  measurements  were  able  to  be  obtained.  Each  component  of  the 
unique  PIV  system  is  described  in  detail  along  with  the  experimental  setup.  In  addition, 
a  computational  procedure  developed  for  high  dynamic  range  (HDR)  analysis  is  presented 
with  accompanying,  sample  results.  Estimations  of  the  measurement  errors  associated  with 
these  results  are  aiso  given.  Finally,  steps  for  improving  the  quality  of  the  experimental 
data  as  well  as  the  analysis  procedure  are  offered  as  suggestions  for  future  investigations. 


I.  Introduction 

The  work  presented  in  this  paper  is  part  of  a  larger,  collaborative  effort  to  investigate  the  turbulence 
associated  with  jet  noise  generation.  In  addition  to  particle  image  veiocimetry  (PIV),  synchronous  near¬ 
field  and  far- field  acoustic  measurements  were  obtained  to  better  quantify  the  noise  generated  by  supersonic 
flows  emanating  from  converging-diverging  (C-D)  nozzles.  Such  noise  is  associated  with  typical  variable- 
area  nozzles  found  on  modern,  high-performance,  military  aircrafts  and  consequently  is  of  interest  due  to 
concerns  over  noise- induced  hearing  loss  as  well  as  degraded  operational  awareness.  The  cumulative  data 
from  these  studies  provides  temporally  resolved,  synchronous  characterization  of  the  near-nozzle  velocity 
field,  the  hydrodynamic  pressure  field,  and  the  acoustic  far  field.  By  studying  the  noise-generating  features 
of  high- temperature,  shock-containing  jets  using  several  measurement  techniques,  the  hope  is  that  a  better 
understanding  of  near-nozzle  flow  conditions  and  their  impacts  on  jet  noise  radiation  can  be  obtained. 

The  following  sections  provide  a  comprehensive  overview  of  the  experimental  research  that  was  performed 
for  this  paper,  that  is,  the  acquisition  of  time-resolved  (TR)  PIV  data  for  the  near-nozzle  flow  field  encom¬ 
passing  the  collapse  of  the  jet  potential  core  (descriptions  and  preliminary  results  for  the  near-field  and 
far- field  acoustic  investigations  arc  available  in  Murray  ct  al.1).  The  aeroacoustic  motivation  for  this  work 
is  provided  by  introducing  the  problem  of  jet  noise.  Following  such  background  information,  accompanying 
sections  are  given  that  depict  the  experimental  setup  and  the  corresponding,  operational  parameters.  The 
specialized  facility  where  measurements  were  performed  is  described  in  conjunction  with  the  unique  PIV 
system  that  was  used.  It  should  be  noted  that  results  obtained  for  this  experimental  investigation  will  be 
presented  at  a  later  time.  The  aim  of  this  paper  is  to  discuss  the  development  of  a  high  dynamic  range  (HDR) 
processing  scheme  that  will  supplement  conventional  forms  of  PIV  data  analysis.  The  proposed  HDR  algo¬ 
rithm  and  its  validation  using  synthetic  particle  images  are  presented  in  section  IV.  Steps  for  improving  this 
HDR  procedure  are  offered  in  the  final  section  along  with  other  concluding  remarks. 
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II.  Experimental  overview 

Despite  over  60  years  of  research  in  jet  aeroacoustics,  a  limited  understanding  persists  in  regards  to 
turbulent  jets  and  the  mechanisms  responsible  for  turbulent  noise  generation.  Not  only  has  a  universal 
definition  of  turbulent  noise  sources  eluded  theorists,  but  what  constitutes  a  source  mechanism  and  how 
such  mechanisms  can  be  rendered  less  efficient  (as  it  pertains  to  sound  generation)  remain  unknown.  From 
an  experimental  standpoint  attempts  to  model,  measure,  and  control  turbulent  jets  have  been  thwarted  by 
instrumental  constraints.  More  often  than  not,  such  constraints  have  resulted  from  inadequate  sensitivities 
and  Insufficient  frequencies  to  accurately  capture  or  noticeably  Influence  phenomena  of  interest.  Until  a 
better  understanding  of  source  mechanisms  is  achieved,  efforts  designed  to  eliminate  or  even  minimize  jet 
noise  radiation  will  continue  to  be  at  the  forefront  of  aeroacoustic  research. 

For  the  purposes  of  this  paper  in  connection  with  the  Jet  Noise  Reduction  (JNR)  program  of  the  Office 
of  Naval  Research  (ONR),  jet  noise  corresponds  to  the  high -amplitude  sound  generated  by  air-breathing 
propulsion  systems,  namely  low-bypass  turbine  engines.  At  present  it  represents  one  of  the  most  acute  noise 
sources  for  the  Department  of  Navy  and  has  been  linked  to  adverse  biological,  mechanical,  and  environmental 
effects.  Such  effects  Include  but  are  not  limited  to  the  noise-induced  hearing  loss  of  Navy  personnel,  the 
structural  degradation  of  Naval  airframes,  and  the  restriction  of  maintenance,  testing,  and  training  sched¬ 
ules  due  to  noise  pollution  in  communities  surrounding  Naval  installations.2  To  counteract  these  issues  and 
combat  the  problem  of  jet  noise,  attempts  are  being  made  to  realize  and  establish  jet  noise  reduction  tech¬ 
nologies  through  coordinated  science  efforts.  The  multitude  of  experiments  undertaken  in  this  investigation 
represents  oue  such  effort.  As  mentioned,  while  the  overall  project  goal  is  to  obtain  a  benchmark-quality 
data  set  that  includes  timc-dcpcndcnt,  velocity-field  measurements  along  with  synchronized  ncar-ficld  and 
far- field  pressure  signals,  the.  research  presented  in  this  paper  is  solely  concerned  with  the  characterization 
of  the  near-nozzle  velocity  field. 

II. A.  Anecholc  Jet  Laboratory 

Experiments  were  conducted  in  the  Anechoic  Jet  Laboratory  (AJL)  at  the  Jamie  Whitten  National  Center 
for  Physical  Acoustics  (NCPA)  on  the  campus  of  the  University  of  Mississippi.  The  AJL  is  a  small  facility 
purpose  built  for  the  study  of  high-tcmperature,  supersonic  jet  noise.3  To  overcome  the  shortcomings  of 
previous  facilities,  specifically  the  NASA  Langley  Small  Anecholc  Jet  Facility  (SAJF),  the  AJL  Is  designed  to 
allow  for  aspiration  of  the  test  chamber.  Using  a  10,000  standard  cubic  feet  per  minute  (SCFM)  fan,  ambient 
air  can  be  pulled  through  the  facility  at  speeds  of  approximately  1  ft/s  (measured  without  jet  flow  in  the 
anechoic  section).  Because  of  upstream  and  downstream  stagnation  chambers,  the  air  actually  percolates 
into  the  l9-by-20-by-8  ft  test  chamber  (measured  between  the  wedge  tips)  through  50%  porosity  sliding 
panels  in  the  upstream  and  downstream  walls.  This  mode  of  operation  results  in  a  very  even  temperature 
distribution  throughout  the  room  while  maintaining  an  acoustically  anechoic  environment.  By  aspirating 
the  entire  chamber,  problems  associated  with  localized  heating  can  be  minimized  along  with  adverse  effects 
on  the  jet  entrainment  due  to  the  enclosed  space.  Figure  1  provides  a  view  of  the  test  chamber  in  the  AJL 
with  various  measurement  systems  in  place. 

The  jet  rig  visible  in  figure  1  and  shown  specifically  in  figure  2  is  supplied  air  from  an  1100  hp  Ingersoll- 
Rand  Centac  compressor  through  a  desiccant  dryer  system.  A  maximum  volumetric  flow  rate  of  5000  SCFM 
of  dry  (-40° F)  air  at  125  psia  enables  continuous  operation  of  the  facility  at  desired  test  conditions.  Control 
valves  operated  in  a  closed-ioop  system  allow'  the  exit  Mach  number  to  be  maintained  within  1%  of  a 
specified  value.  Heat  can  be  added  to  the  flow  through  the  use  of  a  gaseous  propane  burner  system  as 
shown  in  figure  2(b).  The  actual  propane  combustor  is  housed  well  upstream  of  the  nozzie  assembly  and  is 
followed  by  a  ceramic  flow  straigliteuer  and  settling  chamber.  Although  multiple  nozzle  assemblies  exist  for 
this  system,  only  the  configuration  showm  in  the  schematic  that  includes  the  ccntcrbody  section  was  utilized 
for  the  work  In  this  paper.  Investigatioas  to  characterize  the  effects  produced  by  the  other  nozzle  assemblies 
are  currently  underway.  The  compilation  of  results  for  each  of  these  configurations  will  be  presented  in  a 
future  work. 
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Figure  1.  Experiment*]  setup  In  the  test  chember  of  the  Anecholc  Jet  Laboratory  (AJL). 


(a)  Jet  rig  in  the  AJL  (b)  Jet  rig  »chematlc  (Murray  et  a!.1) 

Figure  2.  AJL  propene  burner  system  end  nossle  assembly. 
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II. B.  MHz  frame  rate  PFV  system 

A  MHz  frame  rate  PIV  system  was  developed  through  the  combined  use  of  a  pulse  burst  laser  and  a  high¬ 
speed,  gated  intensified  CCD  framing  camera.  For  its  ability  to  acquire  sequences  of  16  images  at  MHz 
frame  rates,  the  system  allowed  temporally  resolved  velocity-field  measurements  to  be  obtained  for  a  high- 
temperature,  supersonic  jet.  Each  component  of  the  unique  PIV  system  is  explained  in  detail  along  with 
the  experimental  setup.  Although  this  system  was  synchronized  with  both  near-field  and  a  far-field  pressure 
measurement  devices,  only  the  PIV  system  is  considered  for  the  purposes  of  this  paper. 

II. B.  1.  Pulse  burst  laser 

As  has  been  described  in  previous  publications,4-6  a  pulse  burst  laser  system  developed  at  Auburn  University 
allows  a  specified  number  of  high-energy,  MHz  rate  laser  pulses  to  be  formed  for  a  given  burst  of  low- 
energy,  short-duration  pulses.  It  should  be  noted  that  several  upgrades  have  been  made  to  this  system  since 
these  publications,  including  a  new  JDSU  NPRO  126  continuous-wave  (CW)  Nd:YAG  laser  to  enhance  the 
pulse-to-pulse  stability  of  each  burst.  This  component  in  particular  results  in  more  consistent  illuminations 
between  images  and  thus  better  results  in  the  PIV  cross-correlations.  In  addition  to  the  CW  laser,  three 
supplementary  amplification  stages  have  been  Incorporated  Into  the  system  (for  a  total  of  six  amplification 
stages)  to  increase  the  overall  energy  available  for  each  burst  and  consequently  each  individual  laser  pulse. 
A  schematic  of  the  upgraded  pulse  burst  laser  system  is  shown  in  figure  3. 


Pulse  Generation  Pulse  Energy  Amplification  Freq.  Conversion 


Wave  plate 
Optical  Isolator 


Biconvex  lens 
fcassJ  Plano-concave  lens 


- .  Mirror 

0  Polarizer 


Beam  dump 


Optical  crystal 


Figure  3.  Puts*  burst  laser  system  utilised  for  t lm*-resolved  (TR)  PIV  (top  view). 


The  design  of  the  pulse  burst  laser  system  can  be  divided  into  three  fundamental  parts  as  indicated  in  the 
schematic:  the  pulse  generation,  the.  pulse  energy  amplification,  and  the.  frequency  conversion.  As  the  name 
suggests,  the  function  of  the  pulse  generation  stage  is  to  slice  the  output  of  the  CW  laser  into  a  burst  of 
low-energy,  short-duration  pulses.  Slicing  is  achieved  through  the  use  of  an  acousto- optic  modulator  (AOM) 
that  relies  on  the  principles  of  the  acousto-optic  (AO)  effect.  In  particular,  a  piezoelectric  transducer  is  used 
to  produce  acoustic  waves  inside  an  optical  crystal  such  that  the  traveling  waves  cause  variations  in  the  index 
of  refraction  of  the  crystal.  To  an  optical  beam,  these  variations  appear  as  a  sinusoidal  grating  in  which 
the  wavelength  is  equal  to  the  acoustic  wavelength.  By  controlling  when  and  how  frequently  acoustic  waves 
are  produced  inside  the  crystal,  the  generation  of  a  specified  number  of  short-duration  puises  is  possible 
depending  on  how  often  the  CW  input  beam  is  disturbed.  As  with  most  AO  devices,  the  AOM  operates  in 
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the  Bragg  regime  where  most  of  the  incident  light  can  be  diffracted  into  the  first-order  beam  fairly  efficiently. 
This  diffracted  beam  constitutes  the  desired  burst  of  pulses  utilized  in  experimental  applications. 

Following  the  formation  of  low-energy  (nano joule  order),  short-duration  pulses,  the  remaining  stages 
of  the  pulse  burst  laser  consist  of  pulse  energy  amplification  and  frequency  conversion.  Amplification  is 
provided  by  six  flash  lamp-pumped  Nd:YAG  rod  amplifiers  of  increasing  diameter  and  is  necessary  if  the 
pulses  are  to  be  used  for  flu  id- mechanical  measurements.  The  first  three  amplifiers  are  used  in  a  double-pass 
arrangement,  whereas  the  final  three  allow  only  for  single  passes.  Without  going  into  detail,  wave  plates 
and  polarizers  provide  the  necessary  means  for  achieving  double  passes  through  the  first  three  amplifiers. 
Optical  Isolators  between  each  of  the  first  five  amplification  stages  prevent  problems  associated  with  parasitic 
lasiug  and  amplified  spontaneous  emission  (ASE).  By  the  end  of  the  amplification  chain,  pulse  energies  have 
increased  by  a  factor  of  more  than  lO^  lO8  and  generally  reach  leveis  in  excess  of  50  mJ /pulse*. 

The  final  stage  of  the  pulse  burst  laser  system  is  the  conversion  of  the  beam’s  wavelength  from  1064  nm 
to  532  nm.  This  conversion  is  achieved  via  a  nonlinear  process  inside  a  KTP  crystal  and  results  in  an 
unavoidable  loss  of  pulse  energy.  Nevertheless  the  beam,  now  in  the  visible  spectrum,  can  be  used  for 
fluid-meclianical  measurements  including  PIV  and  other  flow  visualization  applications. 

II.B.2.  Cordin  222-  JG  high-speed  camera 

Images  are  acquired  using  a  Cordin  222-4G  gated  intensified  CCD  framing  camera  that  is  capable  of  recording 
16  images  at  a  maximum,  equally  spaced  rate  of  2,500,000  frames  per  secondb.  Such  images  are  captured  with 
a  2048  x  2048  px2  resolution,  although  the  true  resolution  is  less  due  to  the  inherent  intensification  process. 
The  camera  is  able  to  achieve  extremely  high  acquisition  rates  because  it  contains  eight  independently 
controlled  optical  pathways,  each  incorporating  a  microchannel  plate  (MCP)  for  signal  intensification  and 
ultimately  terminating  with  a  Kodak  KAI-4022  CCD  sensor.  Schematics  of  the  camera,  including  an  interior 
view  that  illustrates  four  of  the  optical  pathways,  are  shown  in  figure  4.  By  allowing  each  CCD  to  record 
2  images,  16  total  images  can  be  acquired  over  a  user-specified  time  period.  Furthermore,  because  each 
pathway  is  independently  operated,  temporal  spacing  between  frames  is  variable  and  can  be  set  in  an 
asynchronous  fashion.  Such  flexibility  even  allows  eight  simultaneous  exposures  to  be  made.  This  feature 
is  desirable  since  it  enables  eight  theoretically  identical  particle  images  to  be  obtained,  with  any  differences 
being  directly  attributable  to  error.  More  discussion  on  this  topic  will  be  given  in  the  follow-up  paper 
that  Includes  the  experimental  results.  For  this  work  it  is  sufficient  to  note  that  because  the  camera  can 
acquire  16  images  over  a  user-specified,  extremely  short  time  period,  temporal  resolution  Is  possible  for  all 
captured  fluid  motions.  Additionally,  the  ability  to  obtain  several  particle  images  at  varying  time  intervals 
relative  to  one  another  has  provided  the  means  of  performing  HDR  PIV.  Such  measurements  offer  significant 
improvements  over  conventional  PIV  results  since  optimum  temporal  separations  can  be  selected  for  different 
particle  locations  depending  on  the  local  velocity. 


(a)  Exterior  view  (courtesy  of  Cordin  Company,  Inc.)  (b)  Interior  view  (courtesy  of  Cordin  Company,  Inc.) 

Figure  4.  Side  profiles  of  the  Cordin  222-4G  camere  utilised  In  the  TR  PIV  Investigation. 


“This  value  is  measured  after  the  frequency  conversion  stage  and  thus  accounts  for  the  loss  in  energy  associated  with  doubling 
the  frequency  of  the  Nd:YAG  laser  beam  via  a  KTP  crystal. 

bThis  rate  assume*  a  necessary  CCD  transfer  time  of  3.2  ps  (specified  by  Cordin)  to  ensure  that  the  second  exposure  does 
not  include  images  from  the  first  exposure. 
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II.C.  Experimental  arrangement 

The  combined  use  of  the  pulse  burst  laser  and  the  Cor  din  high-speed  camera  allowed  TR  PIV  measurements 
to  be  made  on  a  high- temperature  (1350°F),  supersonic  jet.  For  all  cases  only  a  smooth  bore  C-D  nozzle 
with  an  upstream  ccntcrbody  section  (shown  in  figures  5(a)  and  5(b))  was  considered  at  over-expanded 
conditions  (Mach  1.55).  Such  conditions  are  typical  of  aircraft  exhaust  during  takeoff  and  low-altitude, 
operation.  The  actual  nozzle  consisted  of  two  conic  sections,  one  contracting  and  the  other  expanding, 
joined  together  to  form  a  supersonic  nozzle  with  a  very  sharp  radius  of  curvature  at  the  throat.  This  near 
discontinuity  at  the  throat  is  significant  since  it  allows  shocks  to  exist  even  when  the  nozzle  is  operated  at 
fully  expanded  conditions  (Mach  1.74).  Tb  illustrate  this  effect,  a  mean  profile  of  the  near-nozzle  velocity 
field  determined  by  CFD  is  shown  in  figure  5(c).  As  indicated  in  the  profile,  the  streamlined  centerbody 
section  was  positioned  well  upstream  of  the  nozzle  contraction.  It  should  be  noted  that  this  nozzle  assembly 
without  the  centerbody  piece  represents  a  1  /  10th  scale  model  of  the  military  power  setting  for  the  General 
Electric  F414  engine. 


(c)  Mean  velocity  field  determined  by  CFD  for  the  centerbody  configuration  at  fully  expanded  conditions  (Murray  et  al.1) 


Figure  5.  The  centerbody  aectlon  Included  In  the  nonle  assembly  Is  shown  In  (e)  end  (b).  The  shock  structures  existing 
et  ftilly  expended  conditions  ere  evident  In  (c). 

The  field  of  view  for  this  work  was  chosen  along  the  bottom  shear  layer  of  the  jet  and  was  centered  at 
a  distance  14  inches  (7  jet  diameters)  downstream  of  the  nozzle  exit.  This  distance  was  selected  based  on 
previous  measurements  indicating  the  collapse  of  the  jet  potential  core.  The  region  imaged  was  slightly  less 
than  16  square  inches  and  was  illuminated  by  a  laser  sheet  directed  vertically  upwards  and  spanning  axially 
along  the  centerline  of  the  jet.  This  particular  orientation  was  chosen  for  a  variety  of  reasons  including  both 
the  need  to  minimize  disruptions  in  the  anechoic  environment  as  well  as  to  ensure  the  most  direct  observation 
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of  any  shear  layer  without  passing  the  light  sheet  through  the  jet  prior  to  imaging.  This  last  point  was 
especially  important  to  prevent  problems  associated  with  aero-optical  distortions.  Figure  6  illustrates  the 
position  of  the  imaged  region  relative  to  the  nozzle  exit  and  the  jet  centerline  (drawing  not  to  scale). 


Fl^ur*  0.  Position  of  th*  Imaged  region  In  tha  TR  PIV  Investigation  relative  to  tha  nossl*  exit  plana  and  tha  J«t 
cantarlln*  (aid*  vlaw).  Th*  aquar*  raglon  Indicate*  the  camera’s  field  of  view. 


Figure  7.  Experimental  arrangement  for  th*  TR  PIV  application  (side  view).  The  aquar*  raglon  enclosed  by  tha  dashed 
Una  Indicates  th*  camera’s  field  of  vlaw. 


Particle  seeding  for  light  scattering  was  achieved  using  aluminum  oxide  (AI2O3)  particles  nominally 
0.1  fj.ni  in  diameter.  A  nitrogen- pressurized  reservoir  filled  with  these  particles  was  connected  to  the  burner 
system  upstream  of  the  nozzle  assembly  and  immediately  following  the  propane  combustor  and  diffuser, 
respectively.  Four  seeding  tubes  were  attached  around  the  burner  system  symmetrically  to  provide  a  uniform 
seeding  density  throughout  the  jet.  To  alleviate  particle  clumping,  a  spinning  propeller  inside  the  reservoir 
formed  a  cloud  of  aluminum  oxide  particles  that  was  subsequently  dispersed  into  the  particle  seeding  tubes. 
Each  connection  between  the  reservoir  and  a  tube  was  made  in  the  supersonic  portion  of  a  miniature  de 
Laval  nozzle  located  at  the  entrance  to  each  tube.  This  arrangement  ensured  that  any  surviving  particle 
clumps  were  sheared  apart  before  being  injected  into  the  burner  system. 

A  schematic  of  the  experimental  setup  is  shown  in  figure  7.  The  532  nm  wavelength  beam  from  the  pulse 
burst  laser  was  passed  into  the  anechoic  chamber  perpendicularly  to  the  jet  axis  and  opposite  the  location  of 
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the  camera.  A  turning  mirror  attached  to  the  stand  for  the  burner  system  allowed  the  beam  to  be  directed 
downstream  of  the  nozzle  exit  before  encountering  a  1000  mm  biconvex  spherical  lens  and  a  second  turning 
mirror  The  beam  was  then  redirected  vertically  upwards  through  a  cylindrical  lens  to  form  the  laser  sheet 
required  for  light  scattering.  Extreme  care  was  taken  to  ensure  that  this  light  sheet  was  oriented  both 
orthogonally  to  the  axis  of  the  camera  lens  as  well  as  to  the  nozzle  exit  plane.  Additionally,  the  placement  of 
the  spherical  lens  allowed  the  thinnest  portion  of  the  light  sheet  to  persist  across  the  camera’s  field  of  view. 

PIV  measurements  were  obtained  by  synchronizing  the  framing  rate  of  the  camera  with  the  pulse- 
generating  rate  of  the  pulse  burst  laser  system.  The  chosen  rate  for  all  cases  was  1  MHz,  meaning  the 
16  images  acquired  by  the  camera  enclosed  a  temporal  window  spanning  15  fi s.  To  achieve  the  most  con¬ 
sistent  pulse- to- pulse  intensity  within  each  burst,  60  laser  pulses  were  generated  for  a  given  burst,  and  the 
most  stable  16  pulses  were  selected  for  synchronization  with  the  16  camera  frames.  The  duration  of  each 
laser  pulse  was  approximately  20  ns  such  that  no  image  streaking  was  observed.  A  Nikon  Nikkor  F-mount 
70-300  mm  objective  zoom  lens  (f/4-5.6G)  was  used  with  the  Cordin  camera  to  acquire  all  image  sequences. 

III.  The  dynamic  velocity  range  in  TR  PIV 

In  addition  to  the  lack  of  three-dimensionality  available  in  PIV,  the  limited  dynamic  velocity  range 
represents  one  of  the  largest  problems  currently  plaguing  the  technique.7  Although  several  methods  have 
been  developed  to  extend  this  range,  flows  containing  a  wide  velocity  distribution  still  present  a  major 
challenge  to  PIV  algorithms.  Because  particle  displacements  in  these  flows  can  differ  drastically  depending 
on  the  local  velocity,  computational  procedures  must  be  designed  to  account  for  temporal  variations  in 
particle-image  patterns.  This  section  describes  a  few  such  approaches  that  aim  to  increase  the  dynamic 
velocity  range  of  TR  PIV  measurements.  The  need  for  temporal  resolution  in  acquired  image  sequences  is 
explained  first,  and  various  strategies  for  evaluating  these  sequences  are  then  presented.  Notable  techniques 
included  in  the  discussion  are  multi-frame  PIV,  adaptive  multi-frame  PIV,  and  adaptive  multi-step  ensemble 
correlation.  A  method  known  as  multiple  pulse  separation  PIV  is  also  mentioned  for  its  relevance  to  improving 
the  dynamic  velocity  range,  however  this  technique  only  provides  time-averaged  results. 

The  current  interest  in  TR  PFV  arises  from  the  need  to  obtain  local  accelerations  and  thus  instantaneous 
pressure  fields  in  unsteady  flow  environments.  Accurate  measurements  of  velocity  time  derivatives  are 
required  in  order  to  evaluate  the  Lagrangian  accelerations  of  fluid  particles  and  subsequently  Integrate  the 
spatial  field  of  the  pressure  gradient.8  For  the  experiments  of  interest  and  other  applications  involving 
aeroacoustics,  the  accuracy  of  such  measurements  is  even  more  important  since  double  temporal  derivatives 
are  required  to  use  the  PIV  data  with  corresponding  aeroacoustic  analogies.9  As  a  result,  reliable  temporal 
derivatives  must  be  achieved  in  PIV  measurements  if  accurate,  subsequent  flow  quantities  are  to  be  obtained. 

One  approach  for  improving  the  accuracy  of  PIV  results  is  to  increase  the  dynamic  velocity  range  of  the 
data.  According  to  Adrian,10  this  range  is  defined  as  the  ratio  of  the  maximum  velocity  to  the  minimum 
resolvable  velocity.  Likewise  it  can  also  be  expressed  in  terms  of  the  maximum  and  minimum  resolvable 
particle  displacements.  Mathematically  these  relationships  can  be  wrritten  as  follows. 

rv\m  —  *hnax  [  ^^mai  / ^®moi  /.\ 

UVK= - =  — 7 - 7771  = -  V*/ 

[  &Ax  l  *  *  \  &Ax 

In  equation  1,  the  dynamic  velocity  range  is  represented  by  DVR,  and  the  maximum  and  minimum  resolvable 
particle  displacements  are  shown  as  and  cr^xi  respectively.  It  should  be  noted  that  only  velocity 

magnitudes  arc  assumed  in  this  definition.  Thus  if  negative  velocities  occur,  then  umox  is  defined  as  the 
larger  of  the  maximum  positive  velocity  or  the  maximum  magnitude  of  the  negative  velocity.  At  present  a 
DVR  of  approximately  200:1  is  the  standard  for  two-dimensional  PIV  measurements.11  In  applications  where 
ranges  begin  to  exceed  this  value,  the  accuracy  of  vectors  in  the  low- velocity  regions  starts  to  deteriorate.  For 
this  reason  data  acquisition  rates  are  generally  chosen  such  that  only  flow  phenomena  of  interest  are  properly 
sampled.  Unfortunately  this  mode  of  imaging  means  that  vectors  calculated  In  other  flow  regimes  have  the 
potential  of  being  highly  inaccurate.  Consequently  high  dynamic  range  techniques  have  been  developed  to 
improve  the  evaluation  of  temporally  resolved  image  sequences. 

Although  increasing  the  time  interval  between  subsequent  frames  is  not  a  preferred  method  for  Improving 
the  DVR,  a  few  studies  have  managed  to  achieve  satisfactory  results  by  applying  this  approach  locally.  The 
problem  with  using  larger  temporal  separations  is  the  reduction  one  experiences  in  the  signal-to-noise  ratio. 
A  greater  pulse  separation  leads  to  increased  losses  between  in-plane  and  out-of-plane  particle  pairs.  To 
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overcome  this  limitation,  procedures  known  as  multi- frame  (MF)  PIV  have  been  developed  that  utilize 
temporally  resolved,  single- frame  image  sequences  to  locally  optimize  particle-image  displacements. 

One  of  the  first  MF  methods  designed  to  improve  the  performance  of  TR  PIV  is  suggested  by  Fincham 
and  Delerce.12  The  approach  revolves  around  a  hierarchical  processing  scheme  that  considers  the  effects  of 
local  fluid  deformation  calculated  during  successive  passes.  Following  an  initial  correlation  of  two  frames 
with  a  temporal  separation  of  St}  displacement  estimates  are  made  for  the  deformation  and  correlation  of 
frames  separated  by  larger  time  intervals  (2 St,  3 St,  etc.).  By  utilizing  larger  time  intervals  in  subsequent 
correlations,  average  pixel  displacements  are  increased  and  thus  the  overall  DVR  enhanced. 

Pereira  et  al.13  introduce  an  alternative  MF  technique  for  dealing  with  temporally  resolved  image  se¬ 
quences.  Termed  adaptive  multi-frame  (AMF)  PIV,  the  technique  aims  at  the  minimization  of  errors  found  in 
low-velocity  regions  by  adjusting  locally  and  dynamically  the  interframe  time  between  particle-image  pairs. 
As  before,  the  algorithm  operates  on  a  local  basis  such  that  a  constant  level  of  accuracy  is  achieved  for  all 
velocity  ranges.  Hain  and  Kahler14  propose  a  further  development  of  this  technique  by  taking  higher-order 
effects  into  account.  For  this  case  images  that  symmetrically  straddle  a  shared,  central  frame  are  considered. 
By  locally  optimizing  the  particle-image  displacements,  an  optimum  temporal  separation  can  be  chosen  for 
each  interrogation  window  that  minimizes  the  relative  measurement  error. 

A  final  method  designed  to  enhance  the  precision  and  robustness  of  TR  PIV  measurements  is  introduced 
by  Sciacchitano  et  al.18  The  technique,  referred  to  as  adaptive  multi-step  ensemble  correlation  (AMEC),  in¬ 
cludes  aspects  of  the  MF  approaches  described  previously  as  well  as  a  method  known  as  correlation  ensemble 
averaging.17  For  a  short  series  of  recordings  separated  by  a  constant  time  interval,  optimum  temporal  sep¬ 
arations  are  locally  evaluated  based  on  error-minimization  criteria.  Correlation  signals  acquired  at  different 
temporal  s pacings  are  linearly  combined  through  the  use  of  homothetic  transformations.  From  comparisons 
with  state-of-the-art  PIV  processing  techniques,  the  AMEC  method  has  proven  to  increase  the  reliability  of 
measured  vectors  and  to  significantly  reduce  both  precision  and  acceleration  errors. 

Until  recently  it  has  not  been  practical  to  use  MF  methods  iu  high-speed  flows  due  to  limitations  im¬ 
posed  by  laser  and  camera  repetition  rates.  To  avoid  the  use  of  excessive  temporal  separations  in  TR  PIV 
applications,  a  multiple  pulse  separation  (MPS)  technique  is  proposed  by  Persoons  and  O’Donovan.15  In  this 
technique  a  series  of  double-frame  images  with  different  pulse  separations  is  recorded  such  that  a  sequence 
with  the  following  temporal  distribution  is  obtained  {[f,  t  +  njr],  [t  +  St,  t  +  St  +  ri2 r],  ...}.  The  inter¬ 
frame  time  (rf<)  remains  constant,  whereas  the  pulse  separation  time  (r)  grows  according  to  a  monotonically 
increasing  multiplier  (nj,  ri2,  etc.).  Once  a  desired  sequence  has  been  acquired,  vector  fields  for  all  pulse 
separation  values  are  evaluated  using  standard  PIV  algorithms.  A  pulse  separation  optimality  criterion  is 
then  applied  locally  to  compute  a  final  displacement  field.  Because  the  results  encompass  multiple  pulse 
separation  values,  the  DVR  is  dramatically  Increased  compared  to  velocity  fields  achieved  by  conventional 
methods.  Despite  this  improvement,  the  MPS  technique  applies  only  to  mean  flow  and  turbulence  quantities 
since  it  is  unabie  to  provide  temporally  resolved  results. 

IV.  Dynamic  evaluation  via  ordinary  least  squares  (DEVOLS) 

As  noted  in  the  previous  section,  flows  containing  a  wide  velocity  distribution  present  a  major  challenge 
to  PIV  algorithms.  The  reason  is  because  particle  motions  In  these  flows  (particle  displacements  in  recorded 
images  of  these  flows)  can  vary  greatly  dependiug  on  the  local  velocity.  Siuce  the  entire  range  of  flow  velocities 
and  thus  particle  motions  cannot  be  adequately  captured  in  a  single  interframe  time,  the  temporal  spacing 
in  conventional  PIV  applications  must  be  chosen  such  that  only  flow  phenomena  of  interest  are  properly 
sampled.  Tb  overcome  this  problem  and  others  related  to  temporal  variations  in  particle-image  patterns, 
HDR  techniques  like  the  ones  described  previously  are  currently  being  developed. 

This  section  presents  a  novel  HDR  approach  designed  specifically  for  the  experiments  of  interest,  that  is, 
the  characterization  of  the  ncar-nozzlc  velocity  field  in  a  supersonic  jet  using  TR  PIV  measurements.  The 
conceptual  idea  for  the  approach  is  based  largely  on  the  MF  method  developed  by  Hain  and  Kahler.14  Using 
the  correlation  results  of  symmetrically  centered  image  pairs  with  increasing  temporal  separations,  a  single 
velocity  field  can  be  constructed  entirely  from  local  evaluations.  This  approach  differs  from  the  previous  ones, 
however,  in  that  individual  vectors  are  determiued  from  the  combined  iufluence  of  measurements  achieved 
at  multiple  interframe  times.  By  considering  the  results  of  several  image  pairs  in  each  local  evaluation, 
significant  improvements  can  be  made  regarding  measurement  accuracy  and  Individual  vector  quality.  The 
following  subsection  discusses  the  proposed  HDR  processing  scheme,  termed  dynamic  evaluation  via  ordinary 
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least  squares  (DEVOLS),  and  its  implementation  into  the  experimental  analysis.  Validation  for  the  procedure 
is  given  in  the  latter  subsection  by  using  synthetically  generated  images  with  known  particle  displacements. 
The  effects  that  particle  density  and  image  noise  have  on  the  algorithm  are  specifically  addressed. 


IV. A.  Proposed  DEVOLS  processing  scheme 

In  the  experiments  of  interest  the  ability  to  obtain  16  particle  images  at  varying  time  intervals  with  respect 
to  one  another  has  provided  the  means  of  performing  HDR  PIV.  Unlike  conventional  PIV  where  only  one 
temporal  spacing  is  available  for  all  velocity  determinations,  the  multiple  combinations  of  image  pairs  in 
this  investigation  enable  a  single  velocity  field  to  be  constructed  from  the  results  of  several  different  local 
evaluations.  In  spatial  regions  where  little  or  no  particle  motions  are  observed  between  subsequent  frames, 
the  results  of  image  pairs  spanning  greater  temporal  distances  are  also  considered.  Thus  it  is  entirely  possible 
for  vectors  in  the  low-vclocity  regions  of  a  flow  field  to  be  determined  using  the  results  of  all  available  image 
pairs.  In  this  way  the  DVR  is  dramatically  improved  because  velocity  ranges  corresponding  to  a  variety  of 
interframe  times  are  properly  and  simultaneously  sampled. 


Figure  S.  Beale  principle  of  high  dynamic  rang*  (HDR)  PIV  In  ralatlon  to  tha  experiments  of  Interest.  A  alngla  valoclty 
fluid  at  t  can  ba  datarmlnad  by  combining  tha  local  a  valuations  of  varloua  Imaga  pairs.  Dapandlng  on  tha  local  valoclty, 
tha  numbar  of  Image  pair*  utilised  for  a  slngla  vector  avaluetlon  is  varlebl*.  The  increased  dynamic  velocity  range 
compared  to  conventlonel  PIV  allows  Improved  measurements  to  be  obtelned. 


A  schematic  illustrating  the  basic  principle  of  HDR  PIV  in  relation  to  the  experimental  investigation  is 
shown  in  figure  8.  As  mentioned,  sequences  of  16  images  were  able  to  be  obtained  in  which  the  St  between 
subsequent  frames  was  specified  to  be  1  /is.  This  particular  interval  was  chosen  such  that  the  maximum 
flow  velocities  and  their  corresponding  particle-image  displacements  would  produce  optimum  results  in  the 
correlation  analyses  between  consecutive  frames.  By  setting  the  data  acquisition  rate  to  properly  sample  the 
maximum  flow  velocity,  the  slower  velocity  ranges  were  inherently  oversampled  and  thus  all  image  pairs  could 
be  used.  Considering  the  temporal  arrangement  of  all  16  frames  in  the  schematic,  the  velocity  field  located 
between  frames  8  and  9  is  determined  by  considering  the  local  results  of  all  8  image  pairs  symmetrically 
straddling  the  point  designated  t.  It  should  be  noted  that  velocity  fields  at  other  points  in  time  could  also  be 
computed,  however  this  particular  position  allows  the  maximum  number  of  image  pairs  to  be  used  with  central 
finite  differencing.  As  is  evident,  for  regions  of  the  velocity  field  containing  the  highest  local  velocities,  that  is, 
the  largest  particle  displacements  between  consecutive  images,  only  the  image  pair  shown  in  red  is  considered 
for  analysis.  Contrarily  for  regions  containing  little  or  no  particle  motions,  multiple  image  pairs  spanning 
larger  temporal  distances  are  considered.  The  image  pair  shown  in  blue  represents  the  case  of  maximum 
temporal  spacing  and  spans  the  entire  sequence  window  of  15  St  or  15  ns.  By  using  the  information  available 
across  all  16  frames  instead  of  only  consecutive  images,  significant  improvements  in  desired  flow  measurements 
can  be  made.  Revisiting  the  mathematical  relation  discussed  previously,  the  reason  is  because  the  DVR  is 
increased  by  a  factor  of  15  compared  to  conventional  PIV  methods.  The  following  equations  demonstrate 
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this  effect  by  combining  the  individual  correlation  results  for  the  minimum  and  maximum  lnterframe  times, 
to:  15 
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Considering  only  the  relations  denoted  by  the  rectangular  boxes  above  (equations  2  and  5),  an  Improved 
DVR  is  obtained. 
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In  the  DEVOLS  processing  scheme  all  eight  symmetrically  centered  image  pairs  are  correlated  in  an 
initial  step.  This  step  represents  the  most  computationally  expensive  and  time-consuming  portion  of  the 
procedure  since  each  correlation  utilizes  multiple  passes  along  with  several  of  the  latest  PIV  processing 
techniques  (window  deformation  and  subpixel  refinement  schemes).  Tb  perform  the  correlation  analyses, 
Dantec  Dynamics  software  (DynamicStudio  v3.31:  Smart  Software  for  Imaging  Solutions)  is  used.  In  order 
for  a  find  velocity  field  to  be  accurately  constructed  from  local  results,  it  Is  imperative  that  the  same  number 
of  interrogation  windows  be  used  in  each  of  the  eight  correlations.  This  requirement  ensures  that  the  vector 
spatial  locations  in  each  displacement  field  are  positioned  and  scaled  in  a  1:1  ratio. 

Following  the  eight  correlation  analyses,  error- minimization  criteria  are  employed  to  select  the  most 
accurate  vectors  determined  at  each  spatial  location.  Such  criteria  are  user-specified  to  allow  for  increased 
processing  flexibility.  To  avoid  errors  stemming  from  particle  accelerations,  displacement  measurements 
are  restricted  by  a  maximum  displacement  limit.  This  upper  bound  ensures  that  high  velocity  regions  are 
only  assessed  by  image  pairs  spanning  the  smallest  temporal  separations.  Contrarily  it  enables  low  velocity 
regions  to  be  evaluated  by  the  maximum  number  of  image  pairs.  In  most  cases  the  highest  accuracies  are 
achieved  when  particle  displacements  are  required  to  satisfy  the  one-quarter  rule.  This  rule  suggests  that 
in-plane  displacements  should  not  exceed  one-quarter  of  the  interrogation  window  size  used  in  the  correlation 
analyses.18  Although  this  rule  is  rendered  obsolete  by  the  use  of  multi-pass/multi-grid  algorithms  in  the 
correlation  analyses  (except  for  the  initial  coarse  grid),  it  provides  a  reasonable  albeit  rudimentary  condition 
for  the  current  processing  scheme.  As  will  be  described  below,  the  algorithm  at  present  requires  valid 
particle-image  displacements  to  increase  in  a  linear  fashion  over  increasing  temporal  spacings.  Modifications 
to  the  program  in  the  near  future  will  remove  this  assumption  by  Inherently  accounting  for  local  acceleration 
effects.  More  on  this  topic  is  given  in  the  concluding  remarks  section. 

In  addition  to  the  particle  displacement  limit,  a  specified  level  of  sensitivity  is  applied  by  the  DEVOLS 
algorithm  when  considering  the  validity  of  vectors.  This  criterion  is  based  on  the  notion  that  for  negligible 
particle  accelerations,  all  image  pairs  should  in  principle  provide  the  same  velocity  measurement  for  a  given 
spatial  location.  Because  particle  motions  under  a  zero-acceleration  condition  remain  constant  with  time,  a 
linear  trend  is  observed  if  their  displacements  are  plotted  over  time.  Considering  only  the  measurements  at  a 
single  spatial  location  satisfying  the  maximum  displacement  restriction,  a  linear  trend  should  be  observed  if 
they  are  plotted  against  their  corresponding  interfrarae  times.  Therefore  using  ordinary  least  squares  (OLS) 
statistics,  a  linear  regression  line  can  be  fit  to  the  data  in  which  the  slope  is  indicative  of  the  local  velocity. 
A  measure  of  how  well  this  regression  line  fits  the  displacement  measurements  is  provided  by  the  coefficient 
of  determination,  denoted  R2.  Depending  on  the  user-specified  sensitivity  level,  l.e.,  the  minimum-allowable 
R2  value,  displacement  measurements  with  the  largest  residuals  are  simply  rejected  until  either  the  R2  value 
exceeds  the  required  tolerance  or  only  a  single,  default  measurement  remains.  The  default  measurement  In 
the  current  algorithm  corresponds  to  the  particle  displacement  determined  by  the  image  pair  spanning  the 
shortest  to.  Thus  the  default  velocity  measurement  for  each  spatial  location  represents  the  conventional  PIV 
measurement  associated  with  the  overall  data  acquisition  rate.  In  the  experiments  of  interest  this  interval 
is  the  1  /xs  separation  time  between  consecutive,  subsequent  frames.  Only  the  default  value  is  used  at  each 
location  unless  supplementary  measurements  provided  by  the  additional  image  pairs  are  deemed  valid  by 
the  user-specified,  error-minimization  criteria. 

Tb  better  explain  the  DEVOLS  processing  scheme,  consider  the  plot  shown  in  figure  9.  In  this  example 
all  eight  displacement  measurements  for  a  given  spatial  location  are  plotted  against  their  corresponding 
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interframe  times.  Because  32  x  32  px2  interrogation  windows  were  specified  during  the  correlation  analyses, 
a  maximum  displacement  limit  was  set  at  12  pixels.  This  value  is  slightly  larger  than  the  one-quarter 
rule  would  suggest.  Based  on  this  restriction  (depicted  by  the  dashed  red  line  in  the  plot),  five  of  the 
eight  measurements  are  immediately  rejected.  As  mentioned,  imposing  a  displacement  restriction  serves  to 
minimize  effects  caused  by  particle  accelerations.  Such  effects  become  increasingly  pronounced  over  excessive 
particle  displacements  since  the  algorithm  currently  assumes  only  linear  trends.  Considering  the  remaining 
three  measurements,  a  linear  regression  line  is  computed  using  OLS  statistics  with  the  added  constraint  that 
the  line  pass  through  the  origin  (since  the  limit  of  Ax  =  0  as  St  -4  0).  Because  the  minimum  R2  value  was 
specified  to  be  0.975,  the  measurement  with  the  largest  residual  (labeled  Outlier  in  the  plot)  is  rejected. 
After  computing  a  second  OLS  regression  line,  the  new  R2  value  is  found  to  exceed  the  specified  tolerance. 
Thus  the  slope  of  this  final  regression  line  represents  the  local  velocity  determined  for  this  particular  spatial 
location.  The  process  is  repeated  for  all  spatial  locations  until  a  final  velocity  field  is  achieved. 


Figure  9.  Graphical  axplanatlon  of  tha  DEVOLS  processing  schama  applied  at  a  single  spatial  location.  Measurements 
are  deemed  Invalid  based  on  a  maximum  displacement  limit  of  12  pixels  end  e  minimum  R3  tolerance  of  0.975  (both  of 
these  criteria  are  specified  by  the  user).  The  slope  of  the  Anal  OLS  regression  line  is  indicative  of  the  local  velocity. 


IV. B.  DEVOLS  validation 

To  validate  the  DEVOLS  processing  scheme  described  previously,  a  time-resolved  sequence  containing  16  syn¬ 
thetically  generated  particle  images  was  considered.  The  particle  density  (p,mg)  in  each  image  was  chosen 
such  that  on  average  each  interrogation  window  would  contain  24  particles.  The  particle  diameters  were 
allowed  to  vary  at  most  by  1  pixel  from  a  nominally  specified  value  of  5  pixels®.  In  addition,  the  particles 
were  allowed  to  exit  the  field  of  view  based  solely  on  their  in-planc  motions.  The  particles  could  also  vary 
in  depth  within  the  light  sheet  (sheet  thickness  for  the  range  [0  1]  was  set  at  0.3),  although  such  positions 
were  fixed  since  no  out-of-plane  motion  was  permitted.  For  this  initial  case  only  zero-noise  conditions  were 
simulated.  Tb  maintain  consistency  with  the  Cordin  CCD  sensors,  the  image  sizes  were  specified  to  be 
2048  x  2048  px2.  It  should  be  noted  that  the  program  used  to  generate  this  sequence  was  a  modified  version 

c Large  particle-image  diameter*  were  specified  to  better  resemble  the  particles  recorded  in  the  experimental  Images.  As 
mentioned,  the  resolution  in  these  images  Is  inherently  reduced  due  to  tbe  intensification  process  associated  witb  the  camera. 
Consequently  the  minimum  resolvable  particle  diameter  associated  with  this  camera  is  slightly  larger  than  would  be  expected 
for  a  non-intenslfied  CCD  camera  with  a  comparable  sensor  size. 
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of  the  synthetic  particle- image  generator  in  PIVlab,  the  time- resolved  digital  particle  image  velocimetry  tool 
for  MATLAB.19 

The  flow  field  simulated  in  the  synthetic  Image  sequence  was  a  Hamel-Oseen  vortex  centered  at  the 
position  (1024.5,  1024.5).  The  maximum  circulation  in  terms  of  maximum  particle  displacement  between 
consecutive,  subsequent  frames  was  limited  to  8  pixels.  This  value  was  chosen  such  that  the  one-quarter  rule 
would  be  satisfied  for  32  x  32  px2  interrogation  windows  over  a  single  interframe  time.  A  50%  overlap  was 
specified  for  these  windows  during  the  correlation  analyses.  A  vortex  was  chosen  to  validate  the  DEVOLS 
processing  scheme  because  it  provides  a  geometrically  simple  case  of  flow  containing  a  wide  velocity  range. 
Depending  on  the  interframe  time,  portions  of  the  flow  field  are  inherently  under-sampled  or  drastically  over¬ 
sampled  when  investigated  by  conventional  PIV.  This  effect  is  evident  In  the  individual  correlation  results 
presented  in  figure  10.  A  small  St  provided  by  the  central  image  pair  can  resolve  the  high-velocity  regions 


AUoMs  Error 


«(p»)  *(vsc) 


(c)  Vector  plot  for  15  St  (d)  Absolute  error  from  the  analytical  result 

Figure  10.  Velocity  fields  determined  by  the  Image  pairs  spanning  St  and  15  St%  respectively.  Every  third  vector  la 
shown  for  clarity  since  each  field  contains  ov  r  16,000  vector*,  As  le  evident  by  the  corresponding  ebsolute  error  plots. 
In  both  caaee  a  Tingle  Interframe  time  le  Inefficient  to  resolve  the  full  DVR  of  the  flow.  The  contour  color  le  Indicative 
of  the  total  deviation  in  pixels  from  the  enelytlcel  solution. 
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near  the  vortex  core  but  in  doing  cannot  adequately  resolve  the  remaining  low- velocity  regions  (figure  10(a)). 
As  a  result,  particle  motions  in  these  regions  approach  the  limit  of  minimum  resolvable  displacement  and 
thus  the  measurement  accuracy  for  individual  vectors  is  of  poor  quality.  Contrarily  a  large  interframe  time 
provided  by  the  image  pair  spanning  the  entire  temporal  window  of  15  St  is  unable  to  resolve  the  core 
(figure  10(c)).  In  this  case,  however,  the  low- velocity  regions  surroundiug  the  core  and  nearing  the  Image 
edges  are  highly  resolved.  The  contour  plot  presented  alongside  each  vector  plot  corresponds  to  the  total 
deviation,  or  absolute  error,  In  pixels  that  exists  for  that  particular  measurement  (<R  and  15  St ,  respectively). 
The  patterns  visible  in  each  plot  illustrate  the  trends  described  previously.  For  both  cases  the  low  DVR  in 
the  measurements  severely  limits  the  viability  of  conventional  PIV  applications.  It  should  be  emphasized 
that  because  this  investigation  represents  the  ideal  case  of  zero  noise,  the  individual  correlations  perform  very 
well.  Still,  the  patterns  visible  in  the  absolute  error  plots  clearly  illustrate  the  effectiveness  of  the  algorithm. 

Using  the  DEVOLS  method  described  previously,  any  number  of  displacement  results  from  the  eight 
correlation  analyses  can  be  considered  at  each  spatial  location  to  obtain  a  more  accurate  measurement  of 
the  velocity  field.  The  fact  that  multiple  interframe  times  are  able  to  be  used  allows  this  measurement  to 
characterize  a  much  higher  DVR  of  the  flow  than  before.  The  plot  shown  in  figure  11(a)  contains  the  final 
HDR  result.  The  error  minimization  criteria  specified  in  the  DEVOLS  algorithm  for  this  case  were  12  pixels 
for  the  displacement  limit  and  0.975  for  the  R2  tolerance.  Unless  specified,  these  values  should  be  assumed 
for  the  HDR  results  contained  in  the  remainder  of  this  paper.  As  before,  the  total  deviation  in  pixels  from 
the  analytical  solution  is  shown  in  figure  11(b). 


Figure  11.  Velocity  Aeld  end  absolute  error  for  ths  HDR  result,  rsspsctlvsly.  The  HDR  result  characterises  the  DVR 
of  the  flow  much  better  then  the  results  shown  In  flgures  10(e)  end  10(c)  since  It  b  not  restricted  to  e  single  lnterfreme 
time.  The  pettern  visible  In  the  contour  plot  cleerly  Ulustretes  this  feet. 


From  the  plots  above,  the  measurement  accuracy  in  the  HDR  result  is  much  better  than  in  the  results 
shown  previously  for  a  single  interframe  time.  The  obvious  reason  is  because  the  DEVOLS  processing 
scheme  can  locally  utilize  the  displacement  results  of  multiple  image  pairs  and  thus  multiple  interframe 
times  to  more  accurately  sample  the  variety  of  flow  regimes.  For  this  flow  field  only  the  image  pair  spanning 
the  shortest  6t  should  ideally  be  used  to  sample  the  high- velocity  regions  near  the  vortex  core.  Contrarily 
multiple  image  pairs  with  successively  increasing  temporal  separations  should  be  used  to  sample  flow  regimes 
located  at  increasing  spatial  distances  from  the  core.  Considering  the  plot  shown  In  figure  12(a),  this  exact 
trend  is  observed.  The  contour  coloring  scheme  is  indicative  of  the  number  of  displacement  measurements 
(Nyee)  utilized  by  the  DEVOLS  algorithm  to  determine  the  final  OLS  regression  line  at  each  vector  location. 
Because  the  minimum  velocities  In  this  flow  field  were  still  significantly  higher  than  the  minimum  resolvable 
limit  (assumed  to  be  0.1  pixels),  only  half  of  the  image  pairs  were  able  to  be  used  at  any  given  location 
based  on  the  specified  settings  (12  pixel  displacement  restriction  and  0.975  R2  tolerance).  Figure  12(b)  also 
illustrates  this  trend  by  plotting  the  velocity  profile  measured  along  the  central,  horizontal  slice  through 
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the  vortex.  As  is  evident,  the  HDR  result  shown  in  red  almost  exactly  matches  the  analytical  solution, 
whereas  the  individual  correlation  results  show  significant  deviations.  The  results  derived  from  the  image 
pairs  spanning  the  largest  temporal  distances  are  especially  ineffective  at  resolving  the  high-velocity  regions 
near  the  core. 


K(|») 

(b)  Velocity  profile  for  the  central,  horizontal  slice 


Figure  13.  Result*  for  tha  vactor  evaluation  flaid  aa  wall  aa  the  cantrai  vaiocity  proflla,  raapactlvaiy.  Tha  trand*  obaarvad 
In  both  caaaa  Indicate  that  only  tha  Image  pair  apannlng  tha  shortest  temporal  dlatanca  la  capable  of  accurately  rasolvlng 
tha  high-valoclty  region*  near  the  vortex  core.  Contrarlly  multiple  Image  pair*  apannlng  increaalng  temporal  distance* 
are  able  to  resolve  the  low-velocity  regions.  Thus  an  Increased  number  of  displacement  results  can  be  utilised  to 
determine  the  velocity  vectors  loceted  at  Increasing  spatial  distances  from  the  core. 


IV.B.l.  Varying  particle  density 

Tb  characterize  the  effect  that  different  particle  densities  have  on  the  DEVOLS  algorithm,  five  temporally 
resolved  image  sequences  were  synthetically  generated  for  the  Hamel- 0 seen  vortex  previously  described. 
All  parameters  (particle  diameter  and  variation,  allowable  depth  within  the  light  sheet,  zero  out-of-plane 
particle  motion,  maximum  circulation  or  particle  displacement  between  subsequent  frames,  and  zero  noise 
conditions)  were  held  constant  between  the  sequences  except  the  total  particle  number.  In  each  case  this 
total  particle  number  was  chosen  such  that  a  desired,  average  number  of  particles  would  be  found  in  the 
interrogation  windows.  Similarly  as  before,  32  x  32  px2  windows  were  specified  in  the  correlation  analyses 
with  a  50%  overlap.  The  particle  densities  considered  in  this  investigation  and  written  in  terms  of  particles 
per  interrogation  window  are  as  follows:  3,  6,  12,  24,  and  48.  The  HDR  results  shown  previously  for  a 
particle  density  of  24  correspond  to  the  same  image  sequence  utilized  in  this  investigation. 

Considering  the  plot  shown  in  figure  13,  the  particle  densities  chosen  for  this  investigation  did  not  appear 
to  affect  the.  DEVOLS  algorithm  in  any  appreciable  way.  A  few  blips  are  noticeable  for  p,mg  =  3,  however 
the  overall  velocity  profile  still  closely  resembles  the  analytical  result.  One  would  expect  the  measurement 
error  to  increase  for  densities  below  3  due  to  the  lack  of  particles  and  consequent  ambiguity  in  the  correlation 
analyses.  One  would  also  expect  the  error  to  increase  above  48  due  to  the  increased  overlapping  of  particle 
images  in  the  recordings.  Quantifying  such  levels  where  the  DEVLOS  processing  scheme  begins  to  experience 
significant  e-rrors  relative,  to  conventional  PIV  algorithms  is  certainly  a  topic  warranting  further  investigation. 
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Figure  13.  HDR  result*  for  the  Hemel-Oseen  vortex  In  which  the  particle  density  was  varied.  Five  temporally  resolved 
Imege  sequences  were  eynthetlcelly  generated  to  conteln  3,  0,  13,  34,  end  48  perticlee  per  Interrogetlon  window.  The 
velocity  profile  shown  corresponds  to  the  centre!,  horisontel  slice  through  the  vortex. 


IV.B.2.  Varying  noise  conditions 

Two  investigations  were  conducted  to  characterize  the  effects  that  different  noise  sources  have  on  the 
DEVOLS  algorithm.  In  both  cases  four  image  sequences  of  increasing  noise  levels  were  considered  against  a 
baseline  sequence  with  zero  noise.  As  before,  a  Hamel-Oseen  vortex  was  simulated  in  the  particle  motions, 
and  all  parameters  were  held  constant  between  sequences  except  for  the  variable  in  question.  In  both  In¬ 
vestigations  the  chosen  particle  density  was  24  (pin,g  =  24),  and  care  was  taken  to  ensure  that  the  initial 
positions  of  all  particles  remained  the  same.  Thus  all  image  sequences  were  identical  in  terms  of  particle 
density,  particle  distribution,  and  particle  motion,  however  they  differed  in  the  prescribed  level  of  noise. 

Salt  and  pepper  noise  The  first  investigation  was  designed  to  simulate  losses-of- pairs  in  the  corre¬ 
lation  analyses.  This  feat  was  accomplished  by  increasing  the  level  of  salt  and  pepper  noise  in  the  image 
sequences.  Such  noise  affects  the  individual  images  within  a  sequence  differently  by  randomly  turning  a 
number  of  pixels  on  to  the  maximum  intensity  value  or  off  to  the  minimum  intensity  value.  The  number  of 
pixels  affected  in  each  image  is  governed  by  a  specified  noise  density  (pnois«)  applied  to  all  of  the  images  for 
a  given  sequence.  Multiplying  this  value  by  the  total  number  of  pixels  in  each  image  provides  the  total  num¬ 
ber  of  pixels  affected.  Noise  densities  of  0,  0.025,  0.05,  0.075,  and  0.1  were  considered  in  this  investigation. 
Sample  particle  images  for  pno j*-  =  0  (baseline)  and  pnoi*c  =  0.075  are  shown  in  figure  14,  respectively. 

The  plot  shown  in  figure  15  contains  the  HDR  results  of  the  five  image  sequences  considered,  namely 
the  baseline,  case,  as  well  as  the  four  cases  of  increasing  noise.  Similarly  to  before.,  this  profile,  corresponds 
to  the  local  velocity  measurements  along  the  central,  horizontal  slice  through  the  vortex.  As  expected,  the 
number  of  deviations  from  the  analytical  solution  appears  to  increase  with  increasing  noise  density,  although 
all  results  resemble  the  analytical  solution  rather  closely.  The  fact  that  large  deviations  appear  at  different 
positions  for  different  image  sequences  Indicates  that  the  salt  and  pepper  noise  applied  to  each  sequence  was 
sufficiently  random.  Thus  the.  losses-of-pairs  in  the.  correlation  analyses  that  resulted  in  the  observed  spikes 
were  clearly  due  to  the  noise  and  not  some  artifact  in  the  underlying  particle-image  distribution. 
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(a)  Salt  and  peppar  noiae:  final  me  —  0  (b)  Salt  and  p«pp«r  noiae:  finolmm  —  0.075 


Figure  14.  Image  regions  containing  different  levels  of  salt  and  pepper  noise.  Each  region  measures  138  x  138  px3  and 
contains  an  tdantlcal  p article  density  and  partlcla  distribution. 


x(px) 


Figure  15.  HDR  results  for  the  HameWOseen  vortex  In  which  the  level  of  salt  and  pepper  noise  was  varied.  Five 
temporally  resolved  Image  aaquances  were  synthetically  generated  to  contain  noise  densities  of  0,  0.035,  0.05,  0.075, 
end  0.1  (where  the  density  multiplied  by  the  number  of  pixels  per  Image  yields  the  total  number  of  pixels  effected). 
The  velocity  profile  shown  correspoods  to  the  central,  horlsootel  slice  through  the  vortex. 
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The  plots  shown  in  figure  16  correspond  to  the  HDR  result  for  pnou«  =  0.075.  Figure  16(a)  shows  the 
entire  velocity  profile,  whereas  figure  16(b)  contains  an  enlarged  view  of  the  region  indicated  by  the  axes. 
Despite  several  wild  deviations  from  the  analytical  solution  by  the  individual  correlation  results,  the  HDR 
curve  follows  the  analytical  result  rather  closely.  As  expected,  the  correlations  of  image  pairs  with  the 
largest  temporal  separations  were  mostly  affected  In  the  high-velocity  regions.  Thus  large  deviations  are 
visible  around  the  peaks  just  outside  the  central  core.  In  contrast  the  image  pair  spanning  a  single  St  was 
mostly  affected  in  the  low- velocity  regions  near  the  image  edges.  Although  the  observed  deviations  for  this 
case  are  significantly  less  than  those  seen  for  the  large  Interframe  times,  they  are  clearly  visible  in  the  zoomed 
view. 


(a)  Velocity  profile  for  the  central,  horizontal  slice 


Flgurs  16.  Result*  achieved  for  the  centrel  velocity  profile  by  tbe  Image  eequenc*  with  />»<><■•  “  0.075.  Plots  ars  shown 
thet  contain  the  entire  velocity  profile  (e)  as  well  as  an  enlarged  view  of  the  region  enclosed  by  the  specified  axes 
values  (b).  The  HDR  result  performed  better  than  the  Individual  correlation  results  for  almost  all  cases. 


Gaussian  white  noise  The  second  investigation  was  designed  to  simulate  intensifier  noise  in  the 
correlation  analyses  by  increasing  the  level  of  Gaussian  white  noise  in  the  image  sequences.  For  clarity  the 
term  Gaussian  refers  to  the  distribution  of  noise  values  in  each  image  (how  frequently  a  particular  value 
appears),  and  white  describes  the  flat  shape  of  the  frequency  spectrum.  As  such  the  noise  applied  to  any 
two  images  (separated  in  time)  for  a  given  sequence  is  statistically  independent  or  uncorrelated.  Analogous 
to  the  noise  density  in  the  previous  investigation,  the  mean  (/i)  and  variance  (<72)  specified  in  the  Gaussian 
distribution  provided  the  means  of  controlling  the  noise  level  generated  for  each  sequence.  In  all  cases  a 
zero-mean  was  specified,  and  variances  of  0,  0.025,  0.05,  0.075,  and  0.1  were  considered.  Sample  particle 
images  for  <r2_0  =  0  (baseline)  and  *2=0  =  o  .075  are  shown  in  figure  17,  respectively. 

Presenting  all  results  in  the  same  order  as  before,  the  plot  shown  in  figure  18  contains  the  HDR  results 
of  the  five  image  sequences  considered,  namely  the  baseline  case  as  well  as  the  four  cases  of  increasing  noise. 
Again  this  profile  corresponds  to  the  local  velocity  measurements  achieved  along  the  central,  horizontal  slice 
through  the  vortex.  Analogous  to  the  trend  observed  for  increasing  particle  density,  the  number  of  deviations 
from  the  analytical  solution  appears  to  increase  for  increasing  values  of  variance.  Nevertheless  all  results 
resemble  the  analytical  solution  to  some  extent. 

The  plots  shown  in  figure  19  correspond  to  the  HDR  result  for  <r2  ^  =  0.075.  Figure  19(a)  shows  the 
entire  velocity  profile,  whereas  figure  19(b)  contains  an  enlarged  view  of  the  region  indicated  by  the  axes. 
Although  the  HDR  curve  follows  the  general  trend  of  the  analytical  solution  rather  closely,  the  deviations  in 
this  result  are  much  more  apparent  than  in  the  analogous  HDR  result  shown  for  the  previous  investigation.  In 
addition  to  the  wild  deviations  experienced  by  the  15  St  and  11  St  correlation  results  near  the  core,  small  blips 
are  evident  in  all  of  the  individual  results  over  much  of  the  profile.  Thus  it  appears  that  the  Gaussian  white 
noise,  at  least  for  the  variances  considered,  had  a  much  more  profound  effect  on  the.  individual  correlation 
results  and  consequently  the  DEVOLS  processing  scheme.  As  is  evident  in  the  enlarged  view,  the  HDR 
result  for  this  case  actually  performed  worse  at  several  spatial  locations  compared  to  many  of  the  individual 
correlations. 
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(a)  Gaussian  white  noise:  =  0  (b)  Gaussian  white  noise:  =  0.075 

Figure  17.  Image  regions  containing  different  levels  of  Geusslen  white  noise.  Each  region  measures  128  X  128  px3  and 
contains  en  Identical  pertlcle  density  and  particle  distribution. 


x(px) 

Figure  18.  HDR  results  for  the  Hemel-Oseen  vortex  In  which  tha  level  of  Geusslen  white  noise  was  varied.  Five 
temporally  resolved  Image  sequences  were  synthetically  generated  to  contain  saro-mean,  Geusslen  white  noise  with 
variances  of  0,  0.026,  0.06,  0.076,  end  0.1  (where  the  meen  end  variance  were  constant  for  ell  cases).  The  velocity 
profile  shown  corresponds  to  the  central,  horlsontel  slice  through  the  vortex. 
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(b)  Enlarged  view  of  the  region  indicated  by  the  axee 


Figure  10.  Reeulte  achieved  for  the  centrel  velocity  profile  by  the  Image  aequence  with  <r*  (0  =  0.075.  Plots  arc  ehown 
that  contain  the  entire  velocity  profile  (e)  ae  well  aa  an  enlarged  view  of  the  region  enclosed  by  the  specified  axes 
values  (b).  Tha  HDR  result  performed  worse  at  savaral  spatial  locations  compared  to  tha  Individual  correlation  results. 


The  reason  the  HDR  result  in  certain  situations  does  not  converge  on  the  analytical  solution,  particularly 
in  the  low-velocity  regions,  is  because  the  default  measurement  is  always  considered  (the  measurement 
provided  by  the  image  pair  spanning  St  cannot  be  rejected  in  the  current  algorithm).  This  unfortunate 
scenario  is  the  result  of  an  inherent  flaw  that  exists  in  the  way  the  current  processing  scheme  is  designed. 
By  not  accounting  for  accelerations  and  thus  assuming  only  linear  displacements,  it  is  possible  in  certain 
cases  for  the  image  pairs  spanning  the  largest  interframe  times  to  provide  a  linear  albeit  incorrect  slope  for 
the  particle-image  displacement  with  time.  Consequently  if  the  OLS  scheme  is  allowed  to  reject  the  default 
measurement  iu  these  cases,  a  regression  line  is  fit  to  incorrect  data.  This  situation  is  avoided  in  the  majority 
of  cases  by  forcing  the  algorithm  to  include  the  default  measurement  in  all  vector  evaluations.  Considering 
figure  19(b),  however,  this  solution  also  results  in  the  algorithm  being  significantly  biased  towards  the  default 
point.  Work  is  currently  being  done  to  develop  an  additional  criterion  that  will  allow  the  default  measurement 
to  be  rejected  or  Ignored  in  the  necessary  cases.  Two  such  methods  are  briefly  discussed  in  the  following 
concluding  remarks. 


V.  Concluding  remarks 

The  development  and  validation  of  a  novel  HDR  processing  scheme  to  supplement  conventional  PIV 
algorithms  has  been  presented.  As  mentioned,  this  approach  is  currently  being  used  to  evaluate  the  TR  PIV 
results  from  an  experimental  investigation  regarding  a  high-temperature,  shock-containing  jet.  The  algo¬ 
rithm,  termed  dynamic  evaluation  via  ordinary  least  squares  (DEVOLS),  offers  substantial  improvements 
over  conventional  PIV  measurements  for  its  ability  to  increase  the  dynamic  velocity  range.  Unique  to 
this  approach  is  an  iterative  validation  scheme  than  enables  multiple  displacement  results  to  be  utilized  in 
the  determination  of  an  individual  velocity  vector.  This  feat  is  accomplished  by  fitting  an  ordinary  least 
squares  (OLS)  regression  line  to  those  displacements  satisfying  a  maximum  displacement  criterion  at  a  given 
spatial  location.  The  slope  of  this  line  is  indicative  of  the  local  velocity.  A  user-specified  tolerance,  that 
is,  a  minimum  allowable  R2  value,  dictates  how  strictly  this  linear  regression  line  must  fit  the  data.  Tb 
validate  the  DEVOLS  algorithm,  a  temporally  resolved  sequence  of  synthetically  generated  particle  images 
was  considered  in  which  the  flow  field  surrounding  a  Hamel-Oseen  vortex  was  simulated.  Effects  due  to 
varying  particle  density  as  well  as  varying  noise  conditions  in  this  flow  environment  were  also  characterized 

In  addition  to  achieving  results  for  the  described  experimental  investigation,  a  couple  of  notable  improve¬ 
ments  to  the  DEVOLS  algorithm  are  planned  for  the  near  future.  The  first  involves  transitioning  the  iterative 
validation  scheme  from  an  OLS  regression  model  to  a  weighted  least  squares  approach.  This  modification  will 
allow  increased  emphasis  to  be  placed  on  measurements  corresponding  to  optimum  particle  displacements. 
Stated  differently,  instead  of  treating  all  valid  displacement  measurements  equally,  an  additional  criterion 
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will  be  used  to  determine  the  optimum  temporal  separation  for  a  given  spatial  location.  Increased  weights 
will  then  be  assigned  to  displacement  measurements  obtained  by  image  pairs  that  span  time  intervals  nearest 
this  optimum  value.  In  addition  to  this  modification,  the  derivation  of  a  least  squares  equation  that  accounts 
for  acceleration  is  also  in  the  works.  This  equation  will  enable  regression  lines  to  account  for  higher-order 
trends  In  the  data  and  not  be  confined  to  linear  fits.  Such  a  scheme  will  also  allow  more  data  points  to  be 
considered  at  each  location  which  will  in  turn  improve  the  overall  measurement  accuracy. 
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